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ABSTRACT 


The  aim  of  this  thesis  is  to  have  a  closer  look  at 
classical  numerical  analysis  as  a  tool  for  solving  the 
elastic  wave  equation  in  its  general  case. 

After  topological  considerations  of  the  parameter 
distributions  at  the  boundary,  and  their  relationship  to 
the  grid  and  the  operator  domains,  we  develop  numerical 
schemes  by  the  direct  method  of  the  numerical  analysis 
(i.e.  the  differential  equation  and  the  external  conditions). 
We  also  examine  the  variational  or  energy  method  where  the 
boundary  conditions  are  naturally  introduced  by  Gauss' 

Theorem  in  the  potential  energy  equation.  For  this 
purpose  we  will  develop  also  a  method  of  discretisation 
of  the  Dirichlet  integral,  and  its  bilinear  form.  Compa¬ 
risons  between  those  schemes  show  the  uniqueness  of  the 
development.  A  finite  difference  example  illustrates  the 
P  -  SV  conversions. 
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CHAPTER  I 


Introduction 

The  properties  of  body  and  surface  waves  have  been 
investigated  by  many  researchers.  More  recently, 
numerical  analysis  has  been  a  powerful  tool.  In  their 
celebrated  paper,  Courant,  Friedrichs  and  Lewy  (1928) 
analysed  the  finite  difference  method  but  it  is  only 
after  numerical  techniques  have  been  used  for  some  time 
for  differential  equations  of  physics  and  engineering 
that  they  have  been  applied  to  hyperbolic  and  transient 
elastic  wave  problems  (Alterman  et  al,  1968).  Alterman 
and  coworkers  studied  the  wave  propagation  in  a  fluid 
sphere  (1968) ,  in  a  layered  elastic  media  (1968) ,  in  an 
elastic  sphere  (1970),  and  in  a  quarter  and  three  quarter 
plane  (19  7  0)  . 

Boore  (1970)  investigated  the  SH  and  Love  wave  pro¬ 
pagations  in  the  transition  region  from  ocean  to  continent. 
Landers  and  Claerbout  (1972)  investigated  the  direct  problem 
of  the  potential  wave  equation  in  elastic  media.  Claerbout 
developed  the  parabolic  approximation  of  the  scalar  wave 
equation  for  direct  and  migration  problems.  Finite  dif¬ 
ference  methods  using  finite  integral  transformations  of 
one  or  more  variables  have  been  widely  used  (Solt^/  1978)  . 

Boore  (1972)  studied  SH  wave  propagation  in  hetero- 
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geneous  media  and  used  successfully  an  equation  with 
variable  parameters.  This  method  had  an  instant  success 
since  cumbersome  boundary  conditions  are  no  longer  needed 
and  arbitrary  boundaries  could  be  considered  without 
adding  any  complexity  to  the  problem.  This  method, 
called  "heterogeneous  media",  was  soon  extended  to  the 
elastic  wave  equation  (Kelly  et  al ,  1976)  .  Finite 

element  methods  (Zienkiewicz ,  1971)  have  been  inves¬ 

tigated  by  Lysemer  and  Drake  (1972)  for  the  direct 
problem  and  by  Marfurt  (1977)  for  the  elastic  wave 
equation  migration  problem. 

Recently  the  necessity  of  obtaining  the  maximum 
information  from  the  data  in  seismic  exploration  has 
given  renewed  interest  in  P-SV  and  SH  waves  and 
consequently  to  the  elastic  wave  equation.  The  first 
aim  of  a  researcher  is  therefore  to  develop  a  more 
efficient  and  optimized  algorithm.  Among  the  methods, 
of  course,  the  "heterogeneous  media"  is  the  most 
attractive.  The  inconsistency  in  the  results  and  the 
disparity  with  the  "homogeneous  media"  method  (Kelly 
et  al,  1976)  has  brought  some  doubt  on  this  technical 
legacy.  Then,  considering  the  "experimental"  aspect 
of  the  so-called  "homogeneous  media"  method  (Chapter  II) 
one  wonders  if  it  is  not  necessary  to  solve  the  problem 
before  computation,  rather  than  to  input  the  equation 
with  external  conditions  in  the  machine  and  carry  on 


endless  trials  and  comparisons. 
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The  objective  is  therefore  to  find  a  scheme  where 
the  boundary  conditions  are  well  defined  in  any  point 
of  a  considered  region  for  any  boundary  shape,  provided 
of  course  that  those  boundaries  pass  by  the  grid  knots. 
For  this  purpose  the  classical  numerical  analysis  offers 
the  choice  between  the  direct  method  and  the  method 
based  on  energy  or  variational  principles  (Chapter  II). 

As  in  any  discrete  problems  we  shall  consider  an 
elementary  domain  and  iterate  the  solution  in  all  the 
considered  regions.  The  direct  approach  will  be 
obtained  through  a  Taylor  series  development,  the 
boundary  conditions  being  introduced  in  the  equations 
(method  developed  in  Chapter  III) . 

It  will  then  be  interesting  to  find  a  solution 
by  the  finite  variational  method,  since  the  boundary 
conditions  are  naturally  included  in  the  integral 
equation  (Chapter  IV).  A  comparison  with  the  direct 
approach  can  give  us  a  new  insight  into  the  problem 
as  well  as  on  the  validity  of  the  solutions  obtained. 

Since  the  source  of  P-SV  conversion  is  an  impor¬ 
tant  question  which  can  be  more  clearly  analysed 
through  the  divergence  and  the  curl  of  the  displacement 
vector  (Chapter  IV)  we  will  investigate  those  equations 
under  two  aspects: 

Compatibility  of  the  equations  with  the  dis¬ 
placement  solution  will  give  a  greater  insight 


■ 

. 

. 
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into  different  aspects  of  the  problem. 
Comparison  with  previous  work. 

The  results  of  this  research  have  been  presented 
by  the  author  at  the  50th  Society  of  Exploration  Geo¬ 
physicists  Convention  (Houston,  1980). 


CHAPTER  II 


NUMERICAL  FORMULATION 

s 

2.1  Equations  of  Motion 

The  equations  of  motion  can  be  deduced  from  the 

principle  of  minimum  potential  energy  in  elasticity. 

Let  u.  be  the  components  of  the  displacement  vector  U 

0.  .  the  stresses  components 

ID 

the  body  forces  components  per  unit  volume 

e  .  .  the  strain  components 

ID 

p  the  density. 

The  deformation  components  are  related  to  the  displacement 
vector  by 


e  .  . 
ID 


,  3u .  9u  . 

Ir — i  +  — ii 

2 l3x .  3x. 1 


(2.1.1) 


The  medium  is  elastic  when  there  exist  a  linear  relationship 
between  stresses  and  strain  (Generalised  Hooke's  law, 

Love ,  1927) 


a  .  . 

ID 


C i j  k£  S  kt 


(  i  ,  j  ,  k  ,  Z 


1,2,3) 


(2.1.2) 


Ci  j  kZ 


being  the  elastic  constants  with 


a  .  .  =  a  .  . 
ID  D 1 
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so  that 


i  jk£  k£i j 

This  follows  from  the  existence  of  an  elastic  potential 


2  W 


Cijk£  Sij  e]aZ 


(2.1.3) 


such  that 


a  .  . 

ID 


9w 


9  e  .  . 
ID 


(2.1.4) 


From  (2.1.2)  and  (2.1.3)  we  obtain 


2  W 


eilQll  +  S22a22  +  S33a33  +  ei2Q12  +  e13Q13  +  623Q23 


...(2.1.5) 


We  shall  represent  the  elastic  potential  by  W(U) .  The 
relation  (2.1.3)  shows  that  W  is  a  quadratic  form  in  the 
strain  components.  In  elasticity  theory  it  is  proved 
that  this  form  is  positive  bounded  below  (Miklin,  1967) . 

If  the  elastic  medium  is  isotropic  then  the  Hooke's 
law  reduces  to  the  Lame  equations 


a..  =  X  0  <5  +ue  .  (2.1.6) 

ID  D  ID 


with  0  =  div  U, 
Lame  constants. 


61  the  Krokener  delta. 
D 


X ,  U  are  the 
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Then,  the  equilibrium  equations  can  be  written  in  the  form 


-  D  .  G  .  .  =  p  . 

D  1 D  i 


(2.1.7) 


where  D . 

3 


9 


9x  . 


+ 


being  the  body  forces  per  unit  volume. 
And  from  (2.1.7)  we  have 


AU  =  F 


(2.1.8) 


where  the  operator  A,  in  the  case  of  an  isotropic  medium,  is 


A  E  (X  +  2y ) V  ( V  .  )  -  PVx ( V  x )  (2.1.9) 


When  the  parameters  are  variables,  (2.1.8)  can  be  written 


PD 


tt 


D.TXD.u.  +  2  PD  .  u  .  ] 
1  3  3  i  i 


-  D , [ PD . u .  +  D.u.]  -  K.  =  0 

3  3  1  l  :  1 


(2.1.10) 


or  in  vectorial  form  (Karal  and  Keller,  1959) 


pii  -  V  [  ( X  +  2  p)  V  .  u  ]  +  Vx[  Vxu  ]  -  K  =  0 


(2.1.11) 


g] 
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Those  equations  are  expressed  in  a  variable  parameter  form, 
since  this  form  is  more  general  and  coincides  with  equation 
(2.1.7)  . 

If  the  propagation  occurs  in  the  (x,z)  plane  the  com- 
ponents  of  u  are  u,\vand  the  P-SV  equations  of  motion  are: 


pD  u  =  D  [  X  [D  u  +  D  W1  +  2  PD  u] 

tt  X  X  Z  X 


+  D  [  P  [  D  w  +  D  u]  ] 

Z  X  z 


PD 


tt 


w 


D  [X  [D  w  +  D  ull  +  2  PD  w] 
z  z  x  z  J 


+ 


D 

x 


[p  [d  u  +  D  w]  ] 
z  x 


(2 . 1. 12) a 


(2 . 1. 12)b 


and  for  the  SH  waves 


pD  V  =  D  (PD  V)  +  D  (PD  V)  (2.1.13) 

tt  xx  z  z 

where  the  vibration  is  perpendicular  to  (x,z). 

The  systems  (2.1.12)  and  (2.1.13)  are  independent  of 
each  other.  p,  X,  P  can  be  constant  in  a  region  (^)  or 
represented  by  a  continuous  differentiable  function.  In 
both  cases  we  have  p,  X,  P,  e  (02  "S  .  Besides  these  systems 
are  submitted  to  auxiliary  conditions  which  are; 

-  The  displacement  and  velocity  fields  at  the  initial 
in  s  tan  t . 

-  The  continuity  of  displacement  and  stresses  at  the 

re.  1 

boundary  (F)  (i.e.  when  p ,  X ,  p  £  CC  )  (fig.  3  ,  page  14) . 

We  are  therefore  led  to  the  following  definitions. 


. 


.  •  . I .S) 
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A  homogeneous  medium  is  a  region  R  (R  c  |Rn )  where  the 
parameters  (or  the  physical  properties  of  the  media)  are 
constants . 

When  the  medium  admits  continuous  derivatives  at  each 
point  P  £  R,  (A,  y  ,  £)  e  CC  ^ ,  then  we  will  consider  the 
medium  as  a  transition  zone.  (Example:  variation  of  velocity 

with  depth,  etc.) . 

A  heterogeneous  medium  is  a  region  R  where  the  physical 
properties  do  not  admit  smooth  variations  in  at  least  one 
of  its  points  P.  (i.e.  3  a  point  P  £  R  such  that  p,  A  or 
y  £  CE  )  .  Consequently  this  point  is  an  element  of  the 
boundary  T . 

Very  often  heterogeneous  media  are  assimilated  to 
transition  zones  by  simple  reason  of  convenience. 


10 


2.2  BASIC  CONSIDERATIONS 

2. 2. a  Physical  laws  governing  the  system 

In  this  section  as  well  as  in  section  2.3,  we  will 
try  to  state  in  the  most  concise  terms  the  broad  lines 
of  any  correct  numerical  solution.  Although  it  may  seem 
too  general,  it  will  constitute  the  base  on  which  our 
problems  will  be  solved. 

Let  us  consider  an  n  dimensional  Euclidean  space  JRn 
and  a  real  time  interval  [0,T]. 

If  $7  is  a  region  of  )R  with  boundary  T ,  we  have 

ft  =  ft  +  r 

U  is  the  displacement  vector  defined  in  the  space  of 
scalar  or  vector  functions,  respectively  defined  in  the 
set  of  points 


ft  x  [o ,T] 


or 


ft  x  [0,T]  and  T  x  [0,T] 
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The  equations  of  motion  (2.1.8)  can  be  written  in  the  more 
concise  form: 

AU(P)  +  F  =  0 
P  e  ft  X  [0  ,T] 

U, 9u/8xi  e  L2 (ft)  (2. 2 .  1) a 


where  A  is  the 


operator 


and  F  =  pD  U  -  K(P)  (2.2.1)b 

K ( P )  representing  the  body  forces  per  unit  volume  at  P. 

The  subsidiary  conditions  being: 

-  The  initial  conditions 


U  ( x  ,  0  )  =  f  ( x  ,  0  ) 

U(x,0)  =  g(x,0)  (2.2.2) 

(The  initial  conditions  being  homogeneous  if  the  "source" 
is  introduced  in  (2.1.1)  as  Body  force  ) . 

-  The  boundary  conditions 
Generally  they  are  classed  as  the  following: 

(a)  The  boundary  F  is  fixed 


U 


r 


o 


(2 . 2  .  3)  a 


■ 


12 


(b)  T  is  free  from  external  forces 


-»• 

a  .  .  «n 

il 


=  0 


(2  .  2  .  3) b 


(c)  A  part  of  the  boundary  is  fixed, 


and  the  other  is  "free". 


In  the  wave  propagation  problem  we  will  be  mainly  concerned 
with  the  conditions  (b)  at  the  surface  and  with  the  continuity 
conditions  at  the  boundary 


U 


=  U 


r 


(displacement  continuity) 


(2.2. 


3) 


-> 

0  .  .  •  n 
il 


r 


-V 

Q  .  .  •  n 

il 


(stress  continuity) 


r 


(  2  .  2 . 4  ) 


where  T  ,  T  represents  the  interior  or  exterior  part  of  the 
boundary . 


If  n  is  the  outward  normal  to  the  boundary  the  stress 


continuity  (2.2.4)b  can  have  the  general  form  (Mikhlin,  1964) 


/  a..  D.  U  cos (n,x, ) 

L  .  il  1  i 


i  >  1 


=  y  a..  D.  U  cos  ( n  ,  x .  ) 

r+  i,j  13  3 


r 


.  .  .  (2.2.5) 


or 


I 


i  /  1 


a  .  .  DU 

13  2 


=  0 


(2.2.6) 


r 


I.  1  " 


' 


13 


DU 


where  D  U  =  -r —  is  the  normal  derivative  to  the  surface  bound- 
->  dn 


n 


ary  T,  the  sign  (+)  being  attributed  to  the  exterior  region 
of  the  boundary  and  the  sign  (-)  to  the  interior.  The 
condition  (2.2.6)  can  be  considered  as  a  generalized  Neumann 
condition . 

2.2.b  The  energy  method 

If  we  construct  the  inner  product 


( AU , U )  =  /  U.AUdft 

ft 


(2.2.7) 


we  have  (Betti  formula) 


(2.2.8) 


( n ) 

where  T  is  the  stress  vector  acting  on  the  surface  of  the 


boundary . 


Since 


/  \ 

J  Ut  '  dT  =  0  (under  the  boundary  conditions  2.2.4) 


r 


then  , 


(AU  f  U)  =  2  /  W(U)dft  >  0 
ft 


(2.2.9) 


<T  .sox  «*.,«;  **  -V-r.  <  *  bnwo* 
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m->.  <ti  m+v,, 


_ 

I 

L 

r  ' 

•  -  - 

P 

D  =  Pi  rp,n  h  h ; 


Fig .  2.1. 


Representation  of  an  elementary  region 
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Thus  the  operator  A  is  positive  definite  ^ .  The  problem 
enumerated  can  be  shown  to  be  equivalent  to  finding  the 
minimum  of  the  potential  energy 

/  (W (U)  -  UF) d ft  (2.2. 10) 

ft 

r  ,  UAU  . 

1  •  e  •  mm  j  (— -  -  UF  )  dft  (2.2.11) 

ft 

which  can  be  expressed  as  minimizing  the  functional 

I  =  j  (AU,U)  -  (U,F)  (2.2.12) 

which  is  the  Dirichelet  principle  (Lions  1970) 
which  states  that 


Inf  I ( v)  =  I (U)  U  e  H ' (ft) 


(2.2.13) 


The  solution  v  e  H ' (ft)  being  unique,  and 

where  H' (ft)  is  the  space  of  finite  energy  (Sobolev  space, 
Sobolev,  1953)  which  is  Hilbert  space  with  norm  equal  to 


(1)  Note:  this  property  is  obvious  in  the  case  of  SH  wave 

2 

where  A  =  -V 


then 


-/  UV*  2  Udft  =  /  (Vu)2dft 
ft  ft 


/ 

r 


o 


and 


(AU, U) 


(  U , AU )  >  0 


' 


BttoiiaatH  rid  pr.  ii  i.-ii:  *  -•  ^  ert  *  afi  *  ;V"  i:W 
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v 

0  n 

2  +  V 

3  v 

2 

H'  (fl) 

2  £ 

L  (Q)  i=l 

-H 

X 

<ro 

L2  (ft)_ 

(2.2.14) 


We  note  that  (2.2.13)  is  the  extension  of  the  Dirichlet 
principle  to  a  more  general  class  of  operators  and  boundary 
conditions . 


2.3  PRINCIPLE  OF  NUMERICAL  SOLUTION 

2.3.  Generalities 

The  general  principle  is  of  course  to  operate  in  a 
finite  dimensional  space. 

Let  us  consider  the  open  region  c  jRn  of  boundary  T 

and  a  basis  W.  ,  .  .  .  ,  W  in  H'  (fi)  such  that 

1  m 

-  W^,  . . . ,  are  linearly  independent  Vm 

-  The  linear  finite  combination 

5.  w.  ,  C.  e  R 


are  dense  in  H'  (tt)  . 

o 


The  problem  has  three  equivalent  solutions  (Lions, 1970) 


(i)  the  usual  formulation  given  by  (2.2.1)  ,  (2.2.4)  . 


(ii)  a  Dirichlet  principle  formulation,  the  solution 


U  being  given  by 
m 


inf  I ( V )  =  1(D)  ,  Um  £  Wm  (2.3.1) 

vew 

m 


which  is  the  Ritz  method. 


\ 
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(iii)  or  by  the  equivalent  form 

(A  U  ,W)  =  (f,W)  V  W  £  W  (2.3.2) 

m  m 

V  P  €  (ft)  =  ft  +  T  c  \Rn 

which  is  the  Galerkin  method. 

We  note  that  the  advantage  of  form  (2.1.1)  and  (2.1.2) 

is  that  the  boundary  conditions  are  included  in  the  equation 

to  solve  (see  note  p.  30)  . 

Then,  the  main  questions  raised  are: 

(i)  the  convergence  and  stability 

(ii)  the  choice  of  the  base  (W  ,  . . . ,  W  ) 

q  m 

The  second  question  being  subject  to  the  matrix 
||  A(W^,W^)||  being  sparse  which  implies  that 

(support  W.}  fl  {support  W,}  =  0 
i  3 

where  {w. ,  . . . ,  W  }  are  functions  characteristic  of  separable 

1  m 

sets  such  that  {w  ,  .  .  .  ,  W  }  are  dense  in  H'  (£2) 

q  m 

and  V  W .  e  H ' (ft)  c  l2 (ft) 

l 

We  will  use  those  properties  in  Chapter  III  to  develop 
finite  difference  schemes  in  which  boundary  conditions  are  implied. 


« 
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2.3.b  Well  posed  problem 

From  the  above  section  we  can  say  the  problem  consists 
in  the  mapping  of  the  initial  data  to  the  space  of  admitted 
solution  U. 

Referring  to  Hadamard  (1953)  the  problem  is  well  posed 
if 

(i)  to  every  set  of  initial  data  corresponds  one  and  only  one 
admitted  solution  U.  (i.e.  A  has  a  unique  inverse  A  ^). 

(ii)  the  solution  U  depends  continuously  on  the  initial  data 
which  implies 

-  if  Uq  is  an  initial  state  which  is  not  element  of  the  analytical 

r>u 

solution,  there  exists  an  element  U  of  the  solution  such  that  lb  -U  I  <£ 

o  o  o 

-  The  stability  (as  a  consequence  of  the  continuity) . 
Besides,  the  geophysical  requirements  are: 

-  To  handle  arbitrary  admitted  data 

-  The  solution  has  to  satisfy  external  constraints, 
i.e.  initial  conditions  and  arbitrary  boundaries. 

-  The  error  truncation  has  to  be  well  defined. 

2.4  FINITE  DIFFERENCE  SCHEME 

2 . 4 . 1 .  a 

The  finite  difference  schemes  are  given  by  the 
discretization  of  the  operators  A  and  Dtt* 


■ 


. 


. 
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We  first  note  that  the  decomposition  of  A  has  the  form 


A  =  l  a  .  .  D.  . 
•  ■  i  1  i  D 


(2.4.1) 


where 


D .  .  =  D  D 

il  x .  x  . 

i  1 


(2.4.2) 


Let  us  consider  in  the  plane  (x,z)  a  grid  of  period 
^lf^2^  anc^  a  point  P(x,z,t)  such  that  at  a  given  instant 
t  =  £At .  We  have  (fig.  2) 


lAt  l 

P(x,z,t)  =  P  (rah . , nh  )  =  P 

1  2  m ,  n 


(2.4.3) 


By  a  simple  Taylor  expansion  we  obtain  the  spatial 

(i.e.  semi  discretization) 


D  U 
xx 


=  K  [U 


-  2U  +  U 


2  m- 1 , n  m,n  m+l,n  12  lx 


1  (h?  D4  U) ] (2 . 4 . 4) 


or 


^  1 
D  U  =  D  U  -  - 

XX  XX 


12h 


4  4 

T  (h^  D  U) 
2  lx 


(2.4.5) 


with 


D  U 
xx  m ,  n 


=  U  -  2U  +  U  _ 

m- 1 , n  m , n  m+1 , n 


(2.4.6) 


or  with  the  notation  in  figure  2.2,  page  24 


D  U  =  — r-  ( U  -  2U  +  U  ) 
xx  0  ,  2  3  0  1 


(2.4.7) 


I- 


• 
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By  the  same  way 


D  U 
zz  m,n 


n-  1 


2  U  + 

m  ,  n 


U  n4.l  1 
m ,  n+ 1 


U 


<\> 

D  U  =  D  U 
z  z  z  z  m  ,  n 


—*  1  2  [h4  D4  U] 

12h22 


(2.4.8) 


with 


% 

D  U 
z  z 


2  U  + 
0 


(2.4.9) 


and 


^  1 

D  U  =  — — —  [u  „  ,  +  U  -  U  -  U  1 

xz  4  h ih  2  m+ 1 , n  + 1  m-l,n-l  m+l,n-l  m-l,n+lJ 


4h 


1  r  2  ,  2  2  2 

nr  [hi  h2  dk  dz 


(2.4.10) 


1  2 


% 

=  D  U 
xz 


4h  h 

1  2 


1*1  h2 


2  2 

D  D  U] 
x  z 


(2.4.11) 


with 


% 

D  U 
xz 


[u6  +  us 


(2.4.12) 


We  note  that  the  second  order  finite  difference  operators 

D  ,  D  are  convolution  operators  of  the  form  (1,  -2,  1) 

XX  z  z  ^ 

which  lead  to  the  well  known  Crank-Nicolson  tridiagonal 
form  (Mitchel,  1978;  Richtmeyer  and  Morton,  1967). 
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By  the  same  way  we  have  for  the  time  variable 


D  U 
tt 


(U£  1  -  2Ul  +  Ul+1 ) 


At 


12  At 


[At4  D4  U]  (2.4.16) 


% 

D  ^  U  =  D  U  - 
1 1  tt 


12At 


[At4  d4  u] 


(2.4.17) 


with 


DttU 


-  2U  +  U 


+  1 


(2.4.18) 


2.4.1.b  Consistancy  of  the  operator  form 

From  (2.4.1)  we  can  write 


r  i  i  2222 

A  =  A  +  T  -Tt-  <h  h  D  D  )  (2.4.19) 

h  h  .  i  i  i  n 

i  3 

2 

since  the  operator  D.  is  positive  definite  (2.2.b) 


~~j —  [h2  h2  D2  D2  U]  >0  (2.4.20) 

hh  i  i  i  -]  — 

i  j 

a. 

Then  A  is  a  lower  bound  of  A  when  h^,  h^  ->  0.  Consequently 
if  A  is  a  closed  positive  definite  bounded  below  operator 
so  is  A.  To  check  the  symmetry  in  £2  it  suffices  to  write 
(2.1.11)  in  the  form 


- 


. 
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'  r  «  'V  O. 

[ ( X  +  2y )  D  +  y  .  .  D  )] 
xx  ij  zz 


(X  +  y)  D 


\(X  +  y )  d 


xz 


xz 


^  't  / 

(X  +  2y)D  +  yD 

zz  K  XX, 


U  (P) 


=  -p  D  U  ( P  ) 


(2.4.21) 


where  P (mn.  ,  mn  _ )  £  c  ^  (fig.  2) 

12  h 

where  is  the  space  spanned  by  P^  ^  and  its  neighborhood. 


We  see  that  the 


% 


from  A  A  conserve  the 


symmetry  property  of  A  V  U_^(P)  ,  P  £  (£3^)  when  P  is  a  point 


of  coordinate  m^h  , 


.,  m.h.  and  £X  is  the  subdomain  defined 
li  h 


by  P  ±  h . . 

i 


(2.1.1)  yields 


A  U  +  DU,  =  Kj" 
h  h  tt  h  h 


(2.4.22) 


we  are  then  led  to  the  following  methods 


(Ahuh+1 


V  + 


n£+l  ou£  n£_1 

uh  "  2uh  +  uh 


At 


U h  =  (K£+1,  U  )  (2.4.23) 


or 


<W  V 


U£+1 

h 


(2.4.24) 


where  and  represent  the  displacement  U(P)  P  -  £2 

h  h 

at  the  instant  £At  and  (£+l)At  or  more  succintly 

P  £  c  iRn  x  [OT) 
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(2. 


The  scheme  (2.4.23)  is  an  implicit  scheme,  the  scheme 
.24)  being  explicit. 

The  implicit  scheme  is  unconditionally  stable  (Mitchell, 


1978)  (see  Section  3-10)  . 


'  -•  m  x  **j1T  • 


2  4 


Ufpn  ,n-l) 


3 


7 


u. 


{m-i;n+0 


4 


hi 

h4 

hi 

h2 

2 


1 

Z 


Hnrwi  n-i ) 

0 


1— >  X 


u 


(T-tin+i ) 


u=um(n  Uf  U2=  Um/n+| 


fi  g  .  2.2 


Lattice  used  in  finite  difference  representation 
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2.5  USUAL  FINITE  DIFFERENCE  APPROACH 

The  usual  finite  difference  methods  are  mostly  evalu¬ 
ated  by  numerical  experiments.  Kelly  et  al.  (1976)  investi¬ 
gated  the  two  main  computational  schemes. 

i)  The  "homogeneous  formulation"  where  the  standard 
boundaries  conditions  have  to  be  satisfied  between  the 
two  media  being  considered. 

ii)  The  "heterogeneous  media"  formulation  where  the 
equation  with  variable  parameters  is  used  and  where  the 
only  boundary  prescribed  are  at  the  surface. 

2. 5. a  The  "homogeneous  approach" 


Equations  (2.1.12)  are  written  in  the  following 


form  : 


O/ 

D 


U=aD  U+$D 


2  „  2  't 

U  +  (a  -  3  ) D  W 


zx 


tt 


XX 


z  z 


'U 

D  W 
tt 


z  z 


„  2^  ,2 
w  +  3  d  w  +  (a 

XX 


(2.5.1) 


The  B.C.  are  for  a  horizontal  interface 


+ 


+ 


W 


W 


u 


u 


(Displacement  continuity) 


' 
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2,+ 

O 

z  z 


2,. 

a 


z  z 


%  + 

o  =  c 

xz  x: 


(Stress  continuity) 


with 


'V 

a 


=  o 


z  z 


0 


xz 


-  0  at  the  surface. 


(2.5.2) 


The  2,  represents  the  finite  difference  expansion,  and 


2—2  ,2  .2  2, 
a  =  a  w  +  (a  -2  3  ) u 
zz  z  x 


2,  2  r2,  2,  , 
a  =  3  [w  +  u 
xz  x  z 


Then  the  systems  (2.5.1)  ,  (2.5.2)  are  applied  to  each 

,  using  as  expansion  of  media  approach  the 
ous  lines  (fig.  2.3) .  For  details  of  the  method, 
see  Alterman  and  Karal  (1968) . 


2.5.b  The  "heterogeneous  media  approach" 


The  equation  (2.1.12)  are  written  taking  into  account 


that  the 


are  variable ,  \Je  then  obtain 


2,  2  22,  22,  2  2, 

D  U  =  aD  U  +  aDU  +  $D  U  +  3DU  + 
tt  xx  xx  zz  xx 


2  2  2, 

(a  -23  )  D  w 

x  z 


22, 

+  3  D  -W 

Z  X 


(2.5.3) 


. 


2  7 


^  2ro  2%  2^  2°o 

D  w  =  a  D  W  +  D  W+$D  W+8D  w  + 

tt  zz  2z  zz  xzz 


2  2  'V' 

(a  -2R  )  d  u 

3  x 


+  3  D  u 

x  z 


(2.5.3) 


surface  conditions  being  expressed  by 


a  =0 

xz 


a  =0 

xz 


(2.5.4) 


with  the  same  procedure  of  fictitious  line  than  the  above- 
mentioned  . 

We  see  that  the  tentative  solution  consists  in 
ignoring  the  boundary  conditions  by  assimilating  the 


boundaries  to  a  limit  of  transition  zone. 


n,  r  -  I 


m-» 


m 

I 
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n 
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fl-fci 


N  TtR'FA'  C  £  X 

J  f  i  cr*}  ~i  ou~s  2-  /  oe 


(  A ,  /u ;  f  )  i. 


n  -\ 

n-i 

nr5 


Z 


Fig .  2.3. 


The  "homogeneous  media"  method 
2. 3. a  Extension  of  the  media  by 
line  at  the  surface 

The  upper  medium  is  extended  at  the  interface 
Extension  of  the  lower  medium  at  the  interface 


2 . 3  .b 
2 . 3  .  c 
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2.6  PROBLEMS  RELATED  TO  USUAL  FINITE  DIFFERENCE  APPROACH 

Since  the  usual  finite  difference  approach  has  mainly 
empirical  basis,  a  plethora  of  "numerical  experiments" 
have  been  published.  The  main  problems  are  well  defined  by 
Kelly  et  al  (1976) : 

i)  The  "homogeneous  media"  method,  although  taking 
explicitly  into  account  the  boundary  conditions,  gives 
undue  instability  at  the  interface.  Numerous  experiments 
and  tests  have  been  done  on  this  effect  (Ilian,  1975). 

This  is  generally  due  to  a  poor  method  of  satisfying  the 
conditions  at  the  boundary.  Besides  the  method  is  cumber¬ 
some  and  can  only  be  applied  to  simple  cases. 

ii)  The  so  called  "heterogeneous  media"  approach, 
although  more  versatile,  has  the  same  problem  as  the  above 
for  the  free  interface  and  has  poor  record  of  stability  at 
any  boundary.  Besides  these  major  problems  the  question  of 
satisfying  the  equation  boundary  conditions  has  been 
avoided  and  the  approximations  and  errors  are  ignored. 

iii)  The  two  above  methods  give  discrepant  results 
(see  Kelly  et  al,  1976). 

iv)  The  connection  between  those  developments  and  the 
other  numerical  methods  based  on  other  principles  has  not 


an 


been  established. 
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v)  The  treatment  of  the  source  is' cumbersome  and  its 
numerical  evaluation  can  be  simplified. 

Note:  Generally  the  validity  of  variational  principle 

or  energy  methods  are  established  in  the  static  case. 


treat  the  mass  x  acceleration  term  of  the  equations  of 
motion  as  an  additional  body  force  (or  nodal  force) . 

In  this  case,  the  problem  is  equivalent  to  consider¬ 
ing  the  displacement  continuous  and  linear  between  two 
time  samples.  If  F  is  such  a  force  the  work  during  an 
interval  of  time  t-tA/2  and  t+At/2  * 


t+ At/2 


p  _  p 

F  .  (u  2-U  2)  dt  =  F 


mean 


t- At/ 2 


This  work  being  equal  to  the  variation  of  kinetic  energy, 
we  have : 


i  .  e  .  mil  *  (U  2-U 


(AU,  (U^-U  ) 


which  implies 


pu  =  F 


mean 


■ 


CHAPTER  III 


THE  PROPOSED  NUMERICAL  SOLUTION 


The  aim  of  this  chapter  is  to  develop  a  second  order 
finite  difference  scheme  valid  in  anisotropic  media  and 
complicated  structure. 

To  do  so,  we  will  have  first  to  consider  a  basic 
topological  question:  The  parameter  distributions  at 

the  boundary  and  their  specificity  with  regard  to  a 
given  operator. 

3.1. a  Spatial  Discretization  and  Operator  Decomposition 

Let  us  consider  at  a  given  instant  T  =  £At  a  point 
P  (x  ,  x  ,  x  )  e  ft  c  )Rn  with 


x 


X 


X 


1 

2 

3 


mh 

rh 

nh 


1 

2 

3 


An  elementary  region  such  that  P  £  £2  can  be  defined  as 
p  (x  ±  h  ,  x  ±  h  ,  x  ±  h  )  and  since  we  operate  in  a 
discrete,  finite  space  the  point  P  will  have  a  finite 
number  of  points  in  its  neighborhood.  i.e.  if  the  number 
of  spatial  dimensions  is  2  we  have: 


31 


■%  ?  3*20 


■ 
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nbh(P0)=  {Px,  Pg}  C  |D  ( A )  =  ft 


where  ID  ( A )  represents  the  elementary  domain  of  the  operator 
A  (  f ig .  3.1). 

|D  can  also  be  decomposed  into  elemental  subdomains 
ID1  such  that 

2n  * 

ID  =  U  ID1 
1=1 

where  1=1,  II,  III,  IV  if  n  =  2 
1=1,  ...,  VIII  if  n  =  3 

spatial  dimensions 

Since  the  operator  A  can  be  decomposed  in  differential 

•  • 

operators  of  the  form  {a.  .D.D.}1  where  {a.  .j1  is  the 

13  i  3  13 

parameter  corresponding  to  the  elemental  subdomain  p1 
we  can  write 


2 

I 

1=1 


.n 


I  r 

(a.  .D.D.)  U  =  )  ( D . a .  . D . ) 

13  i  3  i  13  3 


U 


Consequently,  in  a  discrete  space  the  "homogeneous  media" 
and  "heterogeneous  media"  formulations  should  be  numeri¬ 
cally  and  algorithmically  identical.  We  will  come  back 
later  to  this  problem. 


' 
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3.1.b  Specific  Distribution  of  the  Parameters 

Let  us  consider  the  application  of  a  differential 

operator^  for  instance^  onto  the  previously  defined  domain. 

If  ]D  ( D  )  is  the  domain  of  the  operator  D  ,  we  can 
x  x 

write : 


D(DIU)  =  {u(P_),  U(P  )}  c  iD1 
x  10 


(3.1.1) 


and  ID  ( D  U)  is  the  restriction  of  the  operator  domain, 
x 

The  property  (3.1.1)  remains  valid  if  we  multiply 


U  by  a  coefficient  X,  then 


x  x 


(3.1.2) 


which  implies  that 


ID  (  X  )  1  =  ID  (  D1  (  u)  )  c  ID1 

X 


(3.1.3) 


Hence,  the  parameter  value  can  be  restricted  to  the 


differential  operator  domain. 


Since  we  have  by  definition 


Measure  (h^) 


U(P1) 


with  h 


1 


(Courant  and  Hilbert,  1961) 


. 
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being  a  Lebesgue  Measure  if  U(PQ)  =  m  =  x  and 
UP1  =  (m  +  l) h  . 


Hence,  we  can  write 


r —  Measure  (P 


(3.1.4) 


x 


Then,  the  relative  change  of  the  grid  size  has  the  same 
effect  as  a  variation  of  parameters. 

These  remarks  are  consistent  with  the  fact  that  in 


2  2  2  2 

the  factors  XAt  /ph  ,  yAt  /ph  are  dimensionless. 


Consequently,  the  values  of  the  parameters  do  not  have  to 
be  unique  in  each  elemental  subspace  ID1,  and  can  have 
different  values  corresponding  to  the  operator  domain 

Hence,  the  boundaries  are  not  limited  to  the  basic 
rectangular  grid  shape  and  can  be  diagonal,  which  allows 
more  precise  contour*. 

So,  eight  parameter  values  can  be  allowed  in  the 
neighborhood  of  P  (fig. 5-'). 


I f  we  call 


xi  =  ai  +  xi)/2 


X  2  “  (X2  +  X2)/2 


X  3  '  (X3  +  X3>/2 


X4  =  ‘X4  +  X4)/2 


(3.1.5) 
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The  coefficients  relevant  to  the  operator  D  ,  for  instance, 

X 

will  be  A  and  A  and  A  .  A  will  concern  the  operator  D  . 

J.  j  z  4  z 

In  the  case  of  mixed  derivatives,  D  D  ,  for  instance,  we 

x  z 

remark  that  the  system  can  always  be  written  in  a  conser¬ 
vation  law  form  (Lax  and  Wendroff,  1960) 


D  (DU)  =  D  (a  )  +  D  (a  ) 
t  t  x  xx  z  xz 

or  D  (DU)  =  D  (a  )  +  D  (Q  )  (3.1.6) 

t  t  Z  ZZ  X  zx 

Consequently,  the  parameters  are  relevant  to  the  domain 
of  the  second  derivative. 

These  remarks  are  particularly  important  since  they 
establish  the  relationship  between 

-  Parameters  and  Grid 

-  Parameters  and  Operator. 

We  note  considering  the  domain  D1  for  instance 


ID 1  (  D  U ) 
xx 

ID 1  (  D  U  ) 
z  z 

ID1  (D  U  ) 
x  z 


u  e  [  0  X)  fl  'D1 


U  £  [  0  z)  fl  ID1 


u  £  ( [ Ox)  X 


[0z)  ) 


,  u 


'  U2 


n  id' 


} 

} 


f 


Ul’  U2’  U6} 


(3.1.7) 


: 


Fig .  3.1. 


Distribution  of  the  parameters  at  the  boundary 
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3.2  TAYLOR  DEVELOPMENT  OF  THE  DERIVATIVES 


Instead  of  developing  the  differential  operators 

independently  of  the  parameters  we  shall  develop  the 

terms  a.  .D.D.U  in  each  subdomain  ID1  c  id. 

13  1  D 

Let  us  suppose  for  convenience  a  =  ^/P  which  is 
the  case  of  the  equation  of  motion  (displacement)  in  a 
liquid . 

We  have,  according  to  fig.  3,2in  the  subspace 
ID1  =  ID ^  at  a  given  instant  t  =  PAt 


1  A  ,  2  ,+  [U1  -  u0]  ,+n 

2  p  PlhlDxxU  =  A1  - tT - AlV 


-1-2  3 

Alhi  3  ,+  hi  4 

6  DxxxU  ^1  24  DxxxxU 


with  ID  (D  U  )  c  ID' 
XX 


(3.2.1) 


I  II 

In  the  elemental  subspace  ID  =  |D  we  can  write 


1  X 

2  p 


P0h  D  u 
2  3  xx 


X 


U 

x 


+ 


A3h3 


3 

D  U 

XXX 


A3h3 

24 


D 


4 

U 

xxxx 


+  .  .  . 


with  ID  (  X/pD  ) 
xx 


C  /D 


II 


(3.2.2) 
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After  obtaining  similar  developments  in  the  subspace  ID 


III 


IV 


and  )D  ,  we  obtain  after  summation  (3.2.2) 


— -  D  U  =  2 

p  XX 


hlPl  +  h3p3  U1[U1  -  U0]/hl  +  X3[u3  '  U0]/h3 


+  [X3V  -  XlDxUl  +  [X3h3  *  XlhllD:  0 


XXX 


-  h  +  x 3h 3 ^ D 


xxxx 


U  +  ... 


(3.2.3) 


with  ID  (  A/pD  )  c  ID 
xx 


and  X,  = 


( A  ^  +  A  ^ )  /2 


pj.  -  (pi  +  pi)/2 


X  3  ■  (X3  +  X3)/2 


(p3  +  P3)/2 


by  the  same  method,  we  obtain  in  |D(A/pD  )  c  )D : 

z  z 


£  DZZU  =  p2h2  l  -p-Jh”'  {X2/h2  tu2  -  U0]  +  A4/h4  [u4  '  U0J 


+  [X4D2U  -  A2Dz„]  +  [X4h2  -  A2h2]Dz2zU 


-  h  fX2hl  + 


XXXX 


U  +  •  •  • 


(3.2.4) 


where  A  „  = 


(A2  +  X2)/2 


P2  =  <P2  +  P2>/2 


X 


(X4  +  X4)/2 


p4  “  (p4  +  P4)/2 
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The  expressions  (3.2.3)  and  (3.2.4)  yield,  At  being  the 
time  increment: 


A  2  2 At  ^ 

P  At  DXXU  ■  iT1-p1  +  h3p3  t[A1/h1[u1  -  u0]  ♦  Vh3[U3  -  “„]] 


A  Z 

+  [A3Dxu  -  A1Dxu]}  +  <A  -  A  ) 0 (h  D . u) 


XXX 


—  o  (  Ah4D  u  ) 

^  xxxx 


(3.2.5) 


A  2^  2At 

—  At  D  U  =  : - : - 

p  zz  h  p  +  h  p 


{  [Vh2[U2  ‘  U0]  +  A4/h2[u4  -  U01] 


+  [A4D2U  *  A2D2U]}  +  °  f(A4  -  Al>h\22»! 


At  .  . , ,  4  , 

+  - 0  0  ( Ah  D  u  ) 

,  Z  Z  Z  Z  2 


(3.2.6) 


Since  AAt  /p  has  in 


the  dimensions  of  a 


(length)  the 


error  of  the  above  expressions 


will  be  of  order 


0(h3)r 
0  (h4) 


if  P  is  a  boundary  point 

if  P  is  an  inner  point 

(or  if  the  boundary  is  parallel  to  the 
domain  of  differentiation) 


' 
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The  expressions  [A  D  U  -  A  D  u]  and  [  A-  D  U  -  A„D  ul  are 

lx  4  z  2  z 

boundary  values  naturally  introduced  in  the  equations; 
we  will  keep  them  in  their  classical  differential  form, 


and  take  into  account  the  boundary  conditions  only  after 


replacement  of  the  second  order  differential  terms  of  the 
equations  of  motion  by  the  above  developments. 

We  note  that  for  SH  waves  or  a  liquid  medium  the 
abovement ioned  boundary  terms  cancel.  This  shows  that 
the  "heterogeneous  media"  development  given  by  Boore  (1972) 
in  the  SH  wave  case,  is  correct,  but  this  is  not  the  case 
for  the  elastic  motion  as  we  shall  see  later. 


By  Taylor  Series  we  have  for  D  (D  )  c  D 

X  z 


+ 


(3.2.7) 


Replacing  in  the  above  expression  U  ,  U  by  their 


developments  we  obtain: 


. 
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U6  “  U1  +  U2  -  U0  +  hlh2DxzU 


+  6  [3hlh2DxzzU  +  3hJh2DxxzU]  +  °(h4) 


(3.2.8) 


which  leads  to 


At2  X/p  hlP^Dzxu  =  At2  ^  { [ (u6  -  ux>  -  (u2  -  u 0 )  J 


+  0  (h3)  +  0  (h4)  } (U  e  D1) 


(3.2.9) 


2  —  2^1 

At  X/p  h  p  D  u  =  At  —  {[(u.  -  u  )  -  (u.  -  u  J] 
1  1  zx  lb  0  4 


+  0  (h3) r  +  0  (h4)  }  (U  e  DIV) 


(3.2.10) 


\  + 
2  A3 


At  X/p  h3p3Dzxu  =  At  —  {[(u7  -  u  )  -  (u  -  u  )] 


-  0 (h3)  +  0  (h4)  (U  E  D11) 


(3.2.11) 


At2  X/p  hjPgD  U  =  At2  j/  {[(U0  -  U4>  -  (O  -  Dg) ] 


-  o  (h3) r  +  o  (h4)  (u  e  d111) 


(3.2.12) 


Summation  of  equations  (3.2.9)  to  (3.2.12)  yield: 


■ 


. 
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a  2  A 

At  —  D  u  = 

P  zx 


A  ^ 

2(h1P1  +  l^p  fAl/h2^(u6  "  “l1  "  <U2  '  V  1 


+  U~/h4  (0l  -  u5>  -  (UQ  -  u4)] 


+  A3/h2[(u2  -  Dq)  -  <u7  -  u3)] 


+  A3/h4[(u0  -  u4>  -  (u3  -  Ug) ] > 


+  0  <h3)  +  0  (h4) 


(3.2.13) 


2  ( A  -  A  )  2  3 

The  term  0(h)  is  of  the  form  - - -  At  0(h  )  .  This  * 

h 

term  cancels  if  Uq  is  an  inner  point.  In  this  case 

4 

the  development  is  of  order  0 (h  ) . 

By  the  same  way 


A  2  * 

At  —  D  u  = 


At 


p  xz“  2 (h2P2  +  h4P4) 


{A2/hlt<U6  ‘  V  '  <U1  -  V] 


+  Vh3[<U2  -  U7 '  -  <U0  -  U3>] 


+  A4/h3[(u0  -  o3>  -  (u4  -  u8)] 


X>i[(Ul  -  “o'  -  <U5  -  U4>0 


(3.2.14) 


. 


' 
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If  UQ  is  an  inner  point  (3.2.13)  and  (3.2.14)  reduced  to 


.  2  A  A  2  A 

At  — D  U  =  At  —  D  u 
P  xz  p  xz 


with 


A  .  2  A  A  2  A  r.  . 

At  —  D  u  =  At  — — - - — rl(u  +  u  ) 

p  xz  p  (h2  +  h  )  L  u6  u8 


-  ( U  5  +  U  7 )  ] / ( h  x  +  h  )  +  0(h4)  (3.2.15) 


To  underline  the  boundary  conditions  (3.2.13)  ,  (3.2.14) 

can  be  written: 


At2  -  D  U  = 


At ' 


p  zx  2  (h1P1  +  h 3 p 3  ^ 


U1(U6  ‘  »!>/>>; 


+  X1(U1  *  U5>/h4  -  X3(U7  -  U3>/h2 


+  + 


-  VU3  ‘  V/h4  -(AlDz  -  XlDz)U 


+( X^d+  +  Ad  ) u} 
3  z  3  z 


+  ( A  _  "  A  )0(h2D  U)  +  0(h4)  (3.2.16) 

2  1  3  zz 

h 


. 
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At2  -  D  U  = 


At 


{X*  (U, 


p  xz  2  ( h  2  p  2  +  h4p4)  2  6 


+X2(U2  -  U7)/h3  ‘  X4(U4 


-X>5  ‘  U4)/hl  ‘  X2°xU 


u2)/h 

U8)/h 


1 


3 


«  +  +  .  —  — 
+  X  D  U  +  X  ,D 
4  x  4  x 


} 


+  ( X  -  X  )  0  ( h 2  D  U)  +  Oh4  (3.2.17) 

4  2  xx 

the  expressions  (3.2.16)  ,  (3.2.17)  being  equivalent  to 

(3.2. 16)  ,  (3.2.17). 

For  convenience  if  FD  (  . )  represents  the  terms  of  the 
finite  difference  development  independent  of  the  deriva¬ 
tives;  the  equations  (3.2.5)  ,  (3.2.6)  ,  (3.2.13)  and 

(3.2.14)  yield 


X  A  2  .X  A  2  2 

—  At  D  U  =  F . D ( —  At  D  U) 
p  xx  p  xx 


- 2  - -  [X  D  u  “  X  D  +  u] 

h1P1  +  h3P3  3  x  lx 


+  0  ( h  )p  +  0  ( h  ) 


(3.2.18) 


■ 


' 
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X  .  2  A  2 

-  At  D  U  =  F  .  D  ( —  At  D  U) 
P  zz  p  zz 


2  At  r,  +  , 

V2  +  V4  1  A^u] 


+  0  (h3) r  +  0  (h4) 


(3.2.19) 


-  At2D  U  =  F  .  D  (—  At2D  u) 
P  xz  p  xz 


V4  [A^  '  A^]U 


+  0  (h2) r  +  0 (h4) 


(3.2.20) 


X  *  2  .  2 

—  At  D  U  =  F .  D  ( —  At  D  U) 
P  zx  p  zx 


At' 


hlPl  +  h  3  P  3 


tAlDz  -  A3Dz]u 


+  0 (h2)  +  0  (h4) 


since  consequently  to  section  3.1 


A.D  =  (A.D  +  \+D+)/2 
1  z  1  z  1  z 


X  _  D  =  (A  +  D  +  +  A~D  ) /2  etc. 
3  z  3  z  3  z 


(3.2. 21) 


(3.2.22) 
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3.3  FINITE  DIFFERENCE  DEVELOPMENT  OF  THE  EQUATIONS  OF 

MOTION  (DISPLACEMENT  VECTOR)  WITH  ARBITRARY  BOUNDARIES 

From  (2.1.12)  we  have  in  the  (x,z)  plane: 


D 


D 


U 


tt 


1 1 


(X  +  2y)/pD 


(X  +  2y)/pD 


U 

xx 


W 
z  z 


+ 


+ 


X 

P 


D  W 
zx 


+ 


y 

p 


x 

-  D  U  + 
p  xz 


y 

p 


(D  W  +  D  U) 
xz  z  z 


(D  U  +  D  W) 
zx  XX 


(3.3.1) 


(3.3.2) 


The  auxilliary  conditions  being  given  by 


4* 

u  =  u 


(3.3.3) 


[(X  +  2y )  D  U  +  XD  w]+ 

X  z 

[  (X  +  2y )  d  w  +  Xd  u] + 

Z  X 

(yD  W  +  yD  U)+ 
x  z 


[  (X  +  2 y ) D  U  +  XD  W]  (3.  3.4) 
X  z 

(vertical  boundary) 

[(X  +  2y)D  w  +  Xd  u]~  (3.3.5) 

z  x 

(horizontal  boundary) 

(yD  W  +  yD  u) "  (3.3.6) 

x  z 


where  the  sign  +  designs  the  exterior  part  of  the  boundary. 

We  note  that  condition  (3.3.3)  expresses  the  displace¬ 
ment  continuity,  the  normal  and  tangential  stress  continuity 
being  respectively  expressed  by  conditions  (3.3.4)  ,  (3.3.5) 

and  (3.3.6).  We  note: 

Vertical  boundary  conditions  are  expressed  by  (3.3.4),  (3.3.6) 

Horizontal  boundary-  conditions  are  expressed  by  (3.3.5),  (3.3.6). 


' 


The  free  boundary  case  (i.e.  surface  condition)  is 
expressed  by 

a  =  o  ) 

zx 

"  (vertical  free  boundary) 

G  =0 
XX 

O  =0 
z  z 

\  (horizontal  free  boundary) 

(with  G  =  G  ) 

0=0  zx  xz 

xz 


[(X  +  2p )  D  U  +  XD  w]  =  0  (3.3.7) 

X  z 

[(X  +  2y)D  W  +  XD  u]  =  0  (3.3.8) 

z  x 


yDW  +  ]iDU  =  0  (3.3.9) 

x  z 


The  initial  conditions  and  the  source  problem  will  be 
considered  in  the  next  section. 


Replacing  relations  (3.2.18)  to  (3.2.21)  into 
equations  (3.3.1)  and  (3.3.2)  yield  for  the  general  case: 


"N 
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F.D(D  U)  = 

F.d{(A  +  2y)/pD  u  +  —  D  W 

XX  p  zx 

+ 

y  [  D  W  +  D  ul } 
xz  z  z  J 

+ 

Vl  +  h3p3  [<X  +  2y)3DxU  ‘  [A  +  2P]lDxUl 

+ 

2  At  2  r  -  +  , 

v* .  Vi  KV  ■  W> 

+ 

At2 

hlPl  +  h3p3  fX3DzW  -  VzW] 

+ 

At2  r 

h2p2  H-  h4p4  KV  -  W] 

+ 

0 (h2)  +  0  (h4)  (3.3.10) 

F.D(D  W)  = 

f.d{(A  +  2y)/pD  W  +  —  D  u 

zz  p  xz 

+ 

y  [d  U  +  D  w]  } 

ZX  XX 

+ 

h2P2A+  h4p4  [(X  +  2V)4DzW  -  (X  +  2p>2DzWl 

+ 

hlP;A:  h3p3  tp3Dxw  -  pi°xwi 

+ 

2 

- - -  [Ad  u  -  A  d  u] 

h  2  p  2  +  h  4  P  4  4  x  2  x 

+ 

2 

- - At— -  [y_D  u  -  y-D  u] 

h  Pi  +  h^p^  ^  z  1  z 

+ 

0  (h2) r  +  0  (h4)  (3.3.11) 

(3.3.11) 


4  9 


F.D(DV)  =  F . D {—  D  V  +  ^  D  V> 
tt  p  zz  p  zz 


2  At  r  -  +  , 

+  Z — A 7 — r -  [y„D  V  -  y  D  V] 

h  2 p  2  +  h4  P 4  4  z  2  z  J 


2  At 


— -  [p3D-V  -  ylD%] 


hlPl  +  h3P3 


+  0  (h3) r  +  0 (h4) 


(3.3.12) 


Substituting  the  boundary  conditions  (3.3.4)  ,  (3.3.5)  and 

(3.3.6)  into  (3.3.10),  (3.3.11)  and  (3.3.12)  yield 


F . D ( D  U)  =  F.d{(A  +  2y)/pD  u  +  —  D  W 
tt  xx  p  zx 


+  y  [d  +  d  u]  } 

xz  zz 


- - r -  [Ad  W  -  A  _D  w] 

hlPl  +  h 3 P  3  1  2  3  2 


- — -  [y  D  w  -  y  D  wl 

h  2 p  2  +  h4p4  4  X  2  x 


+  0  (h2)  +  0  (h4) 


(3.3.13) 
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F  D  (D  W)  =  F.d{(A  +  2p)/pD  W  +  —  D 

zz  p  xz 


+  y  [d  +  d  w ]  } 

zx  XX  J 


At' 


[Ad  u  -  Ad  u] 


h2P2  +  h4P4  2  x  4  x 


At 


h 3  p  3  +  hiPi 


[ViDzu  -  u3dzu] 


+  0 (h2)  +  0 (h4) 


(3.3.14) 


F  D ( D  V )  =  F  d{-  (D  V  +  D  V) } 

tt  p  ZZ  XX 


+  0 (h3) r  +  0 (h4) 


(3.3.15) 


3.4  EXPLICIT  SCHEME 


After  substituting  the  F  d{*}  values  given  by  relations 
(3.2.18)  to  (3.2.20)  into  equations  (3.3.13)  to  (3.3.15), 
and  replacing  the  first  derivatives  by  their  second  order 
development  we  obtain: 


U 


2U-U  1  + 


2  (X  +2p  )  <X  +2p  ) 

{-^TT^-[ui-uo]  +  -  3 


hlPl+h3P3 


h 


(u3-u0)} 


At 


2 (hlPl+h3P3*  "2 


A1  X1 
{—  [(W6-V  +  (VW3>]+  ir[(W!-W5) 


X3 

+  (w0-w4)]  +  —  [(w3-w7)  +  (wQ-w2)] 


[[W  -W  )  +  (W  -W  )] 

4 


' 
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2  At 


y 


y. 


+  *2p2+V4  %  [U2‘U0]  +  hT  [Vuo5> 


At 


y. 


+  2(h2p2+h4p4)  t(»6-»2)  +  t"i-w0> 3 


+  hj  [<«2-W7)  +  (wo-«3^  +  hj  [<W8-W4> 


+  (W3-W0)]  +  —  [(w4-w5)  +  (w0-w1) ] } 


+  0 (h2) r  +  0  (h4) 


(3.4 


W1  =  2W-W  1  +  2At 


( X  +2y  ) 

—  { — - —  ( w  -W  ) 

h  _  p „ +h  p „  h „  K  2  0 


2K2  4  4 


At 


<V2V 

+  - : -  (W  -W  )  }  +  t-7t- 

h4  40  2  h2P2+h4P4 


X2 

{-  [(og-u2)  +  tu1-o0>l 


+  TT3  I'W  + 


X 

+  iT3  [<u8-u4>  +  ‘W1 


X 


h 


[‘VV  +  (VV] 


pAf  ^  ^*1  ^3 

+  u  ^ -  {—  (w  -W  )  +  —  (W  -W  )} 

hlPl+h3P3  hl  1  °  h3  3  0 


At ' 


2<hlPl+h3P3)  "2 


{TT  [<W  +  ‘VV1 


•  1) 


. 
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V 


where 


+  577  [(urV  +  (W3  +  5T  [(VU3> 


CO. 


+  'W1  +  577  [<U3~V  +  <u0'u2>]} 


+  0 (h2)  +  0  (h4) 


(3.4.2) 


1  „  -1  2  At 

=  2 V- V  + 


-  [ -  [v  -V  1 

h  p  0+h  p  Lh  1  2  0J 

2  2  4  4  2 


2 

+  577  +  [577  (vi-vo> 


“3 

+  577  (V3-V0)] 


+  0  (h3)  r  +  0  (h4) 


(3.4.3) 


U 


U 


U 


U, 


U. 


0 

=  u 

m  ,  n 

G 

11 

1 

U 

m  ,  n 

=  U 

m+  1 ,  n 

U2  = 

U 

m  ,  n+ 1 

=  u  . 

m-  1  ,  n 

II 

•^r 

D 

U  . 

m  ,  n-  1 

Un- 1 , m+ 1 

U6  ' 

U 

m+ 1 , m+ 1 

=  U 

m- 1 , n+ 1 

II 

00 

D 

m- 1 , n- 1 

X .  =  (X+  +  X. )/2 

1  11 


y .  =  (y+  +  y . ) /2 

1  1  1 


i  =  1  ,  2 


. 
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We  remark: 

i)  The  boundary  conditions  have  been  incorporated 
to  the  equations  with  a  well  defined  second  order  accuracy 
at  the  boundary. 

ii)  The  parameter  distributions  are  governed  by  rules 
seen  in  section  3.2  and  lead  to  a  unique  development  for 
a  given  accuracy. 

iii)  The  system  satisfies  arbitrary  boundaries  including 
diagonal  boundaries,  as  far  as  those  boundaries  pass  by 
the  grid  points,  and  also  take  into  consideration  velocity 
and  grid  anisotropy. 

iv)  No  need  for  cumbersome  "fictitious  lines"  such 
as  defined  by  the  "homogeneous  media"  method  (Alterman 
et  al  (1968,  1970);  Boore  (1970);  Kelly  et  al  (1976); 

even  in  the  case  of  free  boundary  conditions. 

3.5  STRESS  FREE  BOUNDARY  CASE 

The  equations  are  governed  by  relations  (3.3.7)  to  (3.3.9). 
To  satisfy  those  conditions  as  well  as  the  equations  of 
motion,  it  is  sufficient  to  set  the  relevant  parameters  to 
zero  in  equations  (3.4.1)  to  (3.4.3).  Although  the  problem 
is  automatically  solved  by  equations  (3.4.1)  to  (3.4.3),  we 
will  still  write  the  developments  relevant  to  different 


cases  . 


’ 


. 
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This  problem  has  been  the  source  of  numerous  papers 
since  it  concerns  Rayleigh  waves.  The  equations  are  less 
complicated,  but  we  are  in  the  case  where  there  is  maximum 
error  at  the  boundary  (consequently  maximum  instability 
for  inadequate  systems). 

Two  cases  have  to  be  considered;  the  general  case 
where  the  medium  is  non-homogeneous  on  r+ ,  and  the  most 
common  case  when  the  medium  is  homogeneous  on  T+. 

3. 5.  a  The  Medium  is  Non-Homogeneous  on  T  + 

3.5. a.l  Horizontal  stress  free  surface 


The  boundary  conditions  are  expressed  by  equations 
(3.3.8)  and  (3.3.9) . 

For  a  point  P  at  the  surface  the  values  of 

m  ,  0 

,,  ,  ,  ,  ,  n  ,  —  III  -IV 

(A,  y,  p)  relevant  to  the  closed  domain  D  U  D 

are 


i .  e  . 


( A  ,  y  ,  p)  4  =  o 


( A  ,  y  ,  p)  1  =  0 

( A  ,  y  ,  p)  3  =  0 


where  (A,  y,  p)3  means  A3  =  0,  y3  =  0  ... 


(3.5.1) 


Equations  (3.4.1) 


to  (3.4.3)  become 


« 


' 
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U1  =  2U-U~ 1  +  2it 


hipi  +  h3p3  UX1  +  2pi)/h1[ui-U0] 


+  (V2V/h3(U3-V}  +  — T5 - + 

[(W6-Wl)  +  (W2-W  )] 


A3 

+  t(w3-w?)  +  (w0-w2)]} 


2  At 

+  -y-  WV 
h2P2 

At  ^  P  2 

+  77  ~  {—  {  [  (W  -w  )  +  (W  -W  ) ] 

2h2^2  h  &  2  10 


y2 

+  h“  ^  (W2_W7)  +  ( “W  3 )  ] 


+  0(h  ) 


(3.5.2) 


1  -  1  2Afc2  r  ^  ^2  +  2p2 } 

w  =  2W-W  +  { — r2 —  (w2-wo>} 

it2  r  ^ 2 


5h^  [(U6-U2>  +  (U1-U0>J 


+  sr  [(VV  +  (VU3)]} 


2  At 


W1  y3 

-  {_  (W  -W  )  +  —  (W  -W  )} 

h  i p i  +  h 3  p  3  h  ^  1  0  h3  3  0 


+ 


A  2  03 

■~7T - r  [(u  -U  )  +  (U  -U  )] 

i  , +h  o  _  )  h0  6  1  2  0 


2(hlPl+h3P3 


[<D  -U  >  +  (UQ-U2 ) ] }  +  0 (n)  (3.5.3) 
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v1  =  2V-V  1  +  ^  (V  -V  )  +  0  (h3) 

h2P2  h2  2  0 


(3.5. 


where 


U,  y  ,  p)  3  =  j  (X,  y  ,  p)  * 


(X,  y ,  p)  1  =  j  (X,  y  ,  p)  * 


since 


X  =  (X^  +  X  ) / 2  etc. 


and 


( X ,  y,  p)*  =  (X,  y ,  p)  +  =  (X,  y ,  p) 


if  the  medium  is  homogeneous  on  T  . 

If  we  set  in  (3.4.1)  to  (3.4.3)  the  conditions  at 
the  surface 


(P ,  X  ,  y  )  4 
(  P  ,  X  ,  y  )  ~ 
(  P  ,  X  ,  y )  ~ 


(P  #  X ,  y )  2 
(p ,  X  ,  y)  * 

(p,  X,  y ) +  (media  extension) 


as  we  11  as 


<V  V  =  <W2'  V 


(image  condition) 


4) 


we  obtain  identical  results  than  (3.5.2)  to  (3.5.4). 
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3. 5. a. 2  Vertical  stress  free  surface 

The  boundary  conditions  are  expressed  by  equations 
(3.3.7)  and  (3.3.9) . 

The  parameter  values  relevant  to  the  closed  domain 
_  X  X  . .  1 1 1 

ID  U  ID  are  nulls. 


+ 


l .  e  . 


(A ,  y  ,  p)  3  =  0 


( A  ,  y  ,  p)  4  =  0 


( A  ,  y  ,  p)  2  =  0 


(3.5. 


Equations  (3.4.1)  to  (3.4.3)  become 


U  =  2U-U 


-1  2  At 


+  <x1+2ui)/hi  [Vuo] 


+ 


At 


2  A 


2hipi 


(r^  [  )  +  ] 

h  _  6  1  z  D 


A 


+  ^  [(W1-w5)  +  <w0-w4)]} 


■)  Af  ^  ^4 


At2  ,V2 


+  T7T - ~TT - r  [  (W  _ -W  )  +  (W  -W  )] 

2  h 2  P  2  +h  4  P4 )  h!  62  10 


y 


[  (W4-W5)  +  (wQ-wl)]  +  0(h  ) 


5) 


(3.5.6) 


. 


:  ,  :• 
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W1  =  2 W- W~ 1  +  2At 


h  p  +h  p 
2  2  4k4 


( X  +2y  ) 

{ - — -  [w2-w  ] 


^A4+2p4^ 

+  - 5T -  [«4-w0]} 


At 


2'h2P2+h4P4) 


{[ 


X' 


[(VV  +  (urV] 


X4 

+  h7  [(U4-U5>  +  <VV]> 


2  At 


2  y 


hlPl  “l 


h7  [wrwo] 


At' 


P 


2h  1 P 1 


{  h7  t(U6-Ul>  + 


P 1  2 

+  h~  £  ( U  1~ U  5 }  +  (U0"U4)^  +  °(h  >  0.5.7) 


V1  =  2 V- V  1  +  2At 


y 


—  f  —  r  v  -v  i 

h  p  +h  p  L  h  L  2  0J 


y4  2  At 2  yi 

+  —  [v  -V  ]  ]  +  2 ■  —  (V  -V  ) 

h4  L  4  0JJ  hlPl  h±  K  1  0* 


(3.5.8) 


with 

(X, 

y , 

ii 

CM 

Q_ 

1  <A' 

y #  p) 

(X, 

y , 

P)4  " 

I  (A' 

y /  p) 

and 

(X, 

y , 

p,x  = 

(X  ,  y 

'  P)2 

+ 

2 

+ 

4 


+ 

4 


if  the  media  is  homogeneous  on  F  . 
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3. 5. a. 3  Quarter  plane 

The  parameter  values  relevant  to  the  closed  domain 

-II  ,,  -III  ..  -IV 

ID  U  ID  U]D  arenull 


i  .  e  . 

(X, 

y , 

P)  1  =  o 

(X, 

y  / 

o 

II 

CN 

CL 

(X, 

y , 

P>3  =  0 

(X, 

y , 

o 

II 

1  <5 

CL 

Equations  (3.4.1)  to  (3.4.3)  become 


U 


2U~U'1  +  (Xi+2yi)/hi  [U1'U0] 


At 


X 


2hlpl 


{  5T  t<w6-wi>  +  (w2-”o,] 


+  EX  {  if  'Vo” 


'22 

At2 


y 


2 

+ 


2h2P2  hl 


[(W  -w  )  +  (W  -w  )]}  +  0(h  ) 

h  6  2  1  U 


.  .  .  (3.5.9) 


.  • ;  i 
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W1  =  2W-W  1  + 


V  =  2 V- V  + 


2  At 2 

<A2+2y2) 

h2P2 

h  2 

At2 

^2 

2h2p2 

2  At 2 

hipi 

P-, 

r  [w  -w 

hi  1 

At2 

+ 

yi 

S  [(u 

2hipi 

2  At 2 

h  2  P  2 

y~ (v  -v 
2  2  1 

(VV 


{5T  [‘VV  +  (u2-u0)]}+  o<h2) 


2  At 2  yi  , 

^  ^  (Vl"VOl 


+  0  (h  )p, 


(3.5 


with 

(X, 

y , 

11 

pH 

Q. 

T  <X' 

y»  p) *  , 

( x , 

y  , 

P)2  = 

i  (X' 

p.  p > 2  • 

and 

(X, 

y , 

P>t  ' 

(X,  y 

'  p>2 

if  the  medium  is  homogeneous  on  T  . 


3. 5. a. 4  Diagonal  stress  free  boundary 


If 


( X  ,  p,  p) 3  =  0 


(X  ,  y ,  p) 4  =  0 


the  equations  (3.4.1)  to  (3.4.3)  become 


.  10) 
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U  =  2U-U  + 


-1  2At2  Al+2y 


hlpl 


1  1  r  1 

T— 


2  X 


+  ^7  ^  [(W6-W1>  +  (W2-W0>] 


+  iT  [(W] 

4 


2 At 2  y2  r 

^  5^  [W 


At 


2  y 


2  h  2  P  2  ”1 


{__  [(w6-w2)  +  (W  -»0)] 


H  2  2 

+  —  (W2-W7) ]  +  0(h  )r 


W  -  2W-W  +  — —  [A2+2y2]  [^-Wq] 

h2  P2 


At2  X2 


2h 2  P 2  hl 


[ (U1-U0)  +  (W] 


A2 

+  —  (u  -u  ) }  + 
h3  7  2 


2  At 


2  P2  [WrW0] 


hlpl 


+ 


At' 


2h2hlPl  1 


fy:  [U  -o  ]  +  V{  [n  -o  ]} 


At 


+ 


2h2hlP;  fp3  [U2-U0]]  +  °(h  > 


(3.5.11) 


V  =  2 V- V  + 


-1  2  At 


2 

hipi 


[y  [u  -u  1 
L  M1  1  o J 


+  y 2  [u2-u0i  +  0(h3) 

h2P2 


(3.5.12) 


. 
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3.5.b  Stress  free  boundaries  when  the  medium  is  homogeneous 
on  r  + 


The  expansion  consists  in  a  simplification 
of  the  previous  one.  We  will  consider  only  the  horizontal 
case  since  it  is  the  most  common  and  it  will  be  possible 
to  compare  the  results  to  our  work. 


With 


h ..  =  h _  ,  h_  =  h  . 

13  2  4 


(A,y,p)  =  (A,y,p)  =  (A,y,p)3  =  (A,y,p) 


The  schemes  (3.5.2),  (3.5.3)  yield: 


U1  =  2U  -  U  1  +  (A  +  2y )  [U1  -  2  U  Q  +  U3] 

hxP 


At ' 


2hlh2  6 


A [ (W.  -  W, )  +  (W.  -  w  ) ] 

fo  1  J  / 


+  V2  [U  -  o0] 

h2p 


At 


[y  [  (W  -  w  )  +  (wx  -  w3) ] 


2hlh2p LK2  L  '"6  ”7 


0  (h2) 


(3.5.13) 


. 
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2W 


2  At 

2 


[X  +  2y ] [w2 


] 


At 


2h  h  p 
1  2 H 


A2[(u6  -U7)  + 


(ui  -  U  )] 


+ 


At 


2hih2P 


p[<u6 


Ul>  +  <U3  -  V] 


+ 


+  w3] 


2 

+  0 (h  )  (3.5.14) 

3.5.c  Comparison  with  "numerical  experiment  methods" 

Ilian  et  al  (  1975)  ,  Ilian  and  Lowenthal  (  1976)  ,  Ilian 
( 1 9 7 8 )  studied  experimentally  the  stability  of  different 
combination  of  schemes  in  simple  cases. 

The  finding  is  that  the  two  classical  approximations 
using  fictitious  lines  are  unstable  for  (3/cx  >  .5. 

The  range  of  stability  corresponding  to  different 
combinations  of  derivatives  lead  to  the  following  table 


(Ilian,  1978) . 
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Centered  approximation;  Range  of  stability  3/ot  >  .3 

One  sided  3/ot  >  .350 

Composed  3/ot  >  .375 

New  composed  3/ot  >  .000 

The  last  results  are  normal  since  the  stability  is  independent 
of  3/ot  ratio. 

We  can  state  that  equations  (3.5.13)  ,  (3.5.14)  corres¬ 

pond  to  the  most  stable  scheme. 


* 
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3.6  COMPARISON  WITH  THE  CASE  WHERE  THE  BOUNDARIES  ARE 
ASSIMILATED  TO  THE  LIMIT  OF  THE  TRANSITION  ZONE 
("HETEROGENEOUS  MEDIA"  FORM) 

In  the  displacement  vector  case,  simple  calculus 
shows  that  the  terms  independent  of  derivatives  given  by  a 
Taylor  development  of  second  derivatives  (equations 
(3.2.3),  (3.2.4),  (3.2.13),  (3.2.14))  are  identical  to 

the  case  of  a  transition  zone  when  its  thickness  tends 

to  zero.  As  in  previous  chapters  if  we  call  these  develop¬ 
ments  F . D  (  . )  ,  equations  (3.3.10),  (3.3.11)  and  (3.3.12) 

can  be  written: 


( X  +  2 y )  D  u  -  ( X  +  2y)  D  u 


3  x  lx 


FD (D  u)  =  FD  (  .  )  + 

tt  • 


+  [y4Dzu  -  y2Dzul/(h2P2  +  h3p3) 


(3.6.1) 


( X  +  2y )  D  W  -  ( X+2y )  D  w 


4  z  2  z 


FD (D  w)  =  FD  (  .  )  + 

L»  C 


+  [y  D  u  -  y  D  u]/(h  p  +  h  p  ) 
3x  lx  11  33 


(3.6.2) 


F  D(D  v )  =  F  D  (  .  )  +  (Vi  D  v  -  u  D  v) 

1 1  4  Z  /  Z 


+  (y  D  v  -  y  D  v) 


3  x  lx 


(3.6.3) 


The  boundary  conditions  (3.3.4)  to  (3.3.6)  yield 


■ 


' 
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( X+2y)  D  u 
3  x 


( X  +  2y ) , D  u 
1  x 


A  D  w  -  X  D  w 
1  z  3  z 


(X+2y)  D  w 
4  z 


( X  +  2y )  D  w  =  X  D  u 
2  z  2  x 


A  D  u 
4  x 


y 


D  w 
3  x 


y  D  w  =  y  D  u 
lx  1  z 


y  d  u 
3  z 


y  D  u 
4  z 


y  D  u  =  y  D  w  -  yD  w 
2  z  2  x  4  x 


(3.6.4) 


(3.6.5) 


(3.6.6) 


Replacing  (3.6.4)  to  (3.6.6)  in  equations  (3.6.1)  , 
(3.6.2)  yield 


F  D  (D  u)  =  F  D  (  .  )  +  2  At 

tt 


(Ad  w  -  A  d  w) 
2  1  z _ 3  z 

hl^l  + 


+  2  At 


(A  d  w  -  A  d  w) 
2  2  x  4  x 


tt 


+  2At 


h2P2  + 

h3P3 

)  +  2  At 2 

(A2DxU 

-  A  d  u) 
4  x 

h2P2 

+  h  p 

4  4 

<piV  - 

P3DzU) 

hipi  + 

h3P3) 

(3.6.7) 


(3.6.8) 


F  D  (D  v)  =  F  j)(  .  ) 


(3.6.9) 


(since  for  SH  waves  y  D  v  -  y  D  v  =  0  and  y  D  v  - 

4  z  2  z  3  x 

y iDxv  = 


0 


Vv(P)  p  e  ti) 
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Equations  (3.6.7)  to  (3.6.9)  show  that 

i)  if  P(x,z)  is  an  inner  point  the  developments  are 
of  course  identical  since  the  boundary  terms  cancel. 

ii)  For  SH  waves  the  two  methods  are  identical  since 
the  differential  terms  at  the  boundary  cancel  (eqn. 

3.6.9).  This  lets  us  forecast  different  results  for  P.SV 
waves'  potential  equations.  Although  they  have  the  same 
form, the  boundary  terms  will  not  vanish. 

iii)  For  the  displacement  vector  equations  the  error 
is  considerable  and  equal  to 

At 2  (A  -  X  )D  w/(h  p  +  h  p  ) 

1  3  z  11  3  3 

+  At2(P2  -  P4)Dx«/<h2P2  +  h4p4) 
for  the  horizontal  component  u. 

and  At2(A2  -  X4>Dxu/(h2P2  +  h4P4) 

+  At2(P1  -  y3)Dzu/(h1P1  +  h  3  P  3 ) 

for  the  vertical  component  w. 

Therefore  we  can  see  that  assimilation  of  the  boundary 
conditions  to  the  limit  of  a  transition  zone,  although  correct 
in  the  case  of  SH  waves  (Boor e  1972  )  ,  is  not  acceptable 


* 

■ 
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in  the  case  of  a  displacement  vector  P-SV./  as  in  Kelly  et  al 
"heterogeneous  media"  equations  (1976)  We  shall  consider 
later  the  case  of  P-SV  potential  equations  since  this 
case  has  been  studied  analytically  by  Gupta  (1966)  . 
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3.7  INITIAL  CONDITIONS  AND  SOURCE  REPRESENTATION 

Alterman  and  Karal  (1968)  solved  the  difficulties 
encountered  in  the  neighborhood  of  the  source  by  sub¬ 
tracting  the  displacement  due  to  the  source  from  the 
total  displacement  field  in  a  rectangular  region  surround¬ 
ing  the  source.  The  direct  source  contributions  being 
given  from  the  known  analytical  solution  for  the  source 
in  an  infinite  region,  and  displacement  continuity 
conditions  being  applied  at  the  boundary  of  the  rectangular 
region.  This  method  is  also  applied  by  Kelly  et  al  (1976) . 

Alterman  and  Aboudi  (1970)  ,  Aboudi  (1971)  give  a 
numerical  treatment  of  seismic  sources  in  elastic 
media  which  are  equivalent  to  body  forces. 

We  remark  that  although  the  problems  can  be  super¬ 
posed,  we  shall  distinguish  two  different  problems. 

3.  7. a  ) 

The  initial  displacement  and  velocity  field  are 
prescribed  in  the  considered  region  ft 

U  (P)  =  g  (x  ,  0) 

U  (P)  =  f  ( x  ,  0  ) 

Vp  e  S) 

U  ,  U  EL2  (ft) 


(3.7.1) 


■ 


. 

■ 


where  is  solution  of  the  homogeneous  equation  of  motion 


at 

In  this  case,  no  special  treatment  is  necessary  (see  Boor.e 
1970) . 

3  .  7  .  b  ) 

The  displacement  is  due  to  a  source  S(x,t)  defined 

in  a  region  ft  c  ft  of  duration  [0,t]  such  that 

s 


s (P)  =  o 

VP  ( x , t )  t  x  [o,T)  (3.7.2) 
s 

which  is  the  common  case  of  seismic  exploration. 

If  S  (x,t)  represents  the  displacement  potential  the 
equations  of  motion  yield 

Au  +  pD  u  +  pD^, [grad  S  (x,t)]  =  0 

1 1  1 1  p 

ft  X  [o  ,  T )  (3.7.3) 

s 

and  the  source  will  act  as  a  body  force  since  it  is  a 
body  force. 

For  S  (x,t)  representing  the  displacement  source  we 


shall  have 


■ 
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AU  +  P°ttU  ‘  pDttSd(X,t)  =  ° 

in  Q  x  [  0  ,  T  )  (3.7.4) 


We  note  that  equations  (3.7.3)  and  (3.7.4)  express  very 
clearly  the  reciprocity  between  source  and  receiver. 

In  the  discrete  domain  (3.7.4)  becomes 


AU  +  pD  U  -  pD  S  ( x  ,  t )  =  0 

1 1  1 1  d 


The  development  of  D  S^(x,t) 


be  ing 


PDttsd(x't) 


(in  0,  x  [  0  /  T  )  ) 
s 


(3.7.5) 


(3.7.6) 


for  simple  reason  of  causality. 

In  this  case,  no  restriction  has  to  be  imposed  on  the 

source  function,  except  that  for  a  numerical  treatment  S, 

a 

has  to  be  band  limited,  that  is  its  amplitude  spectrum 
has  to  be  non  zero  over  a  finite  range  of  the  transform 
variable,  to  avoid  aliasing  errors  (see  next  chapter). 

Simulation  of  forces  such  as  defined  by  Aboudi  (1971) 
Archambeau  (1968)  ,  Burridge  and  Knopoff  (1964)  ,  Burridge, 
Lapwood  and  Knopoff  (1964)  can  be  applied  without  complica 


t i on s  . 
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We  note  that  for  second  order  source  development  (3.' 
wi 1 1  be 


S1  -  s° 

p°ttVx't)  ^  2p[~A~t~"~ 


At 


where  S  (P),  S,(P)  =  0  for  P  $  SI  x  [o,t] 

d  d  s 


.6) 


X 


K  =S(t)D  Sp(x  ) 


f*  g  3.2 


Fig  .  3.2. 


Representation  of  a  source  with  a  Gaussian 
distribution. 


7  4 


fig  3-3-b 


Fig.  3.3.  Forces  representing  a  source 

3. 3.  a  Dilatation  source  (superposition  of 

dipoles) ; 

3.3. b  Shear  source  (superposition  of. couples). 


3.8  GRID  EDGE  BOUNDARIES 


To  limit  the  area  of  computation  it  is  essential 
to  introduce  artificial  boundaries.  Those  constraints 
are  of  Neumann  type  and  the  grid  sides  behave  as  perfect 
reflectors.  The  problem  can,  of  course,  be  overcome 
expensively  by  widening  the  grid  in  such  a  way 
that  the  undesired  reflections  do  not  disturb  the 
required  solutions. 

If  there  is  at  present  no  perfect  way  to  simulate 
an  absorbing  boundary,  a  number  of  methods  have  been 
developed  to  cancel  the  unwanted  reflections  by  intro¬ 
ducing  some  constraints  at  the  boundaries.  Lysemer  and 
Kuhlemeyer  (1969)  achieved  attenuation  at  the  boundaries 
by  applying  viscous  tractions  such  that 


normal  stresses  =  apv  w 

P 


shear  stresses  =  bpv  u 

s 


(3.8. 


where  u,  w  are  as  usual  the  time  derivatives  of  the 
displacement  vector  components,  v  and  v  the  compres- 
sional  and  shear  velocities,  and  a,b  are  dimensionless 
coefficients.  Cartellani  (1974)  used  a  similar  approach 
where  the  viscous  forces  are  represented  by 


' 

. 
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a  =  -V  pU  •'  n 
n  pK 


T  =  -V  pU  x  n 
n  s 


(3.8.2) 


where  G  ^  and  T  ^  are  respectively  normal  and  tangent  stresses; 

~T~  '  _  , 

u  is  the  velocity  vector,  n  and  t  are  the  unit  normal  and 
tangent  vector. 

Those  viscous  dash-spots  are  at  the  base  of  most 

* 

techniques  used  for  absorbing  boundaries.  The  quality  of 
the  attenuation  is  a  function  of  the  frequency  and  angle 
of  incidence  of  the  incident  wave  (Lysmer  and  Kuhlemeyer 
(1969)  ,  Castellani  (1974)  ,  Robinson  (197  7)  )  . 

W.D.  Smith  (1975)  took  advantage  of  the  phase  difference 
between  free  (Neuman)  and  clamped  (Dirichlet)  boundaries 
to  attenuate  the  reflections  by  adding  them.  A  closer  look 
shows  that  the  method  is  insufficient  for  elastic  waves 
since  Rayleigh  waves  do  not  exist  for  the  clamped  boundaries. 

The  proposed  approach  is  to  simulate  increasing 
internal  friction  in  a  buffer  zone  (Voigt  Solid,  Ewing 
et  al,  195 7 ).  Then,  the  stresses  can  have  the  form 


Q 


1 

xx 


a 


XX 


a 


xz 


1 


a 


xz 
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If  we  simplify  by  writing  a  =  b ,  we  obtain 


1 

0 

xx 


( 1  +  a 


|->o 

d  t  XX 


1 

a 

xz 


(1  + 


(3.8.3) 


Then,  the  equations  of  motion  become 


u 


w 


1 1 


1 1 


(1  + 


a  + 

XX 


d  a  1 
z  xz J 


(1  + 


)  [3  a 

z 


+ 

z  z 


3  a  ] 

X  ZXJ 


(3.8.4) 


(3.8.5) 


when  the  equations  (3.8.4),  (3.8.5)  are  a  simplified 

form  of  Voigt's  relation.  For  further  simplification 
we  admit  that: 


a (U1-u‘1) /2At  =  - £  2  U . 


r* w*  r>u 


U  =  (  1  -  £  )  [  D  0  +DQ  ! 

tt  X  XX  z  XZ 


r*u 


w..  =  (1-£)[dq  +  d  a  ] 


tt 


Z  Z  Z  X  zx 


(3.8.6) 


(3.8.7) 


To  realize  this  model  it  suffices  to  multiply  the  equation 

2  2 

by  a  coefficient  equal  at  k  =  (1  -  £  )  <  1  this  coefficient 

decreasing  progressively  to  avoid  reflections  in  the 
buffer  zone.  The  effect  is  the  same  as  increasing 


. 
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the  grid  size  by 
the  velocities, 
sufficient  as  a 


- - — —  >  1  or  decrease  progressively 

( 1  -  e  ) 

Generally  8  to  12  spatial  steps  are 
buffer  zone. 
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Fig.  3.4.  Absorbing  boundary  conditions  (Castellani  method) 

3. 4.  a  Ratio  of  reflected  and  incident  energy 

for  compre ssional  waves  as  functions  of 
angle  of  incidence 

3.4. b  Energy  ratio  for  shear  waves 


Cas>'e^^av"  ^7*4  ^ 
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3.9  ACCURACY 

The  accuracy  is  determined  by  the  sum  of  the  rounding 
error  and  the  truncation  error.  Let  U(P)  be  the  analyti- 
cal  solution  and  U(P)  be  the  numerical  solution  at  a  point 
P  . 


2 

P (x,Z ,t)  c  ft  C  R  x  [o,t] 

3. 9. a  Rounding  error 

The  rounding  error  is  a  consequence  of  the  computa¬ 
tional  procedure.  We  have  at  the  initial  moment  (Smith, 
1978) 

N°(P)  =  U°(P)  -  r°(P)  (3.9.1) 

0 

where  r  is  the  initial  vector  of  rounding  error  and 
U°  the  initial  numerical  value  (U?P  )  c  u  (P)  )  . 

If  we  consider  the  two  time  step  numerical  recursive 
solution  of  the  form 


U 


j 


'b  j  - 1 
AUJ 


(3.9.2) 


then 
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N 


1 


%  0 

AN 


N 


j 


0 

r 


Jj-1  1 

A  r 


and 


U 


j 


0 

r 


,  Jj-1  1 

+  A  r 


+ 


+ 


(3.9.3) 


(3.9.4) 


This  shows  that  ths  rounding  error  vector-  r  propagates 
the  same  way  as  the  numerical  solution  vector  U. 

If  A  ^  are  the  eigenvalues  of  A1  i  =  1(1)  (N  -  1) 
thenequation  (3.9.4)  converges  if 


max  A  .  <  1 . 

i  1 


3.9.b  Truncation  error 

a. 

Let  F  D (U (P) )  =  0 


(3.9.5) 


represent  the  finite  difference  equation  at  P(x,z,t). 
If  U  represents  the  exact  solution  of  the  differential 
equation  then  the  truncation  error  is  defined  as 


£  (P)  =  F . D ( U ( P ) ) 


with  P  c  []  c  (Rn  x  [0,T]  (3.9.6) 


Consequently  £  is  derived  from  the  Taylor  development 


residual  terms. 


■ 
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3 . 9 . b . 1  Case  of  P.SV  waves. From  (3.9.6),  (3.3.10)  and 

? 

(3.3.11)  We  have  that 


£  = 


A  (A,p)0(h2) 


2 

h  p 


+  "4“  (  A  ,  p ) 0 ( h  4 )  +  0  ( At  2 )  (3.9.7) 

h  p 


A  2 

where  — r—  (A,y)0(h  ) 

h  p 


=  0  if  P  is  an  inner  point 


Then 


£  =  0 ( h 2 )  +  0 ( A t 2 )  if  P  c  ft 


3.9.b.2  Case  of  SH  waves 


From  equation  (3.3.12)  we  obtain 


£  =  yT-  0(h3) 

h  p 


U  4  2 

+  0(h)  +0  (At  ) 

h  p 


(3.9.8) 


but 


Ay  3.  1 

2  0 (h  } r  3 1  (yi 

h  p 


y  -j )  — r  d  v 

3  ,2  xxx 

h 


+  77  <^2 


y  .)  -  D  v 

4  2  z  z  z 

h 


(3,9.9) 


since  the  boundary  conditions  can  be  written 


y  d+u 

1  x 


y  d  u 
3  x 


y  d+u 
2  z 


y„D  u 

4  z 


(3.9.10) 
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we  have 


0  (h3) 

h  p 


0  . 


Then 


e  =  0  (h2 , At2) 


since  (3.9.10)  can  be  written 


( 2  n  + 1 ) +  (  2  n  + 1 ) - 

y  D  =  yoD 

lx  3  x 


(2n+l) + 

y  D 


D (2n  +  l)  - 


(3.9.11) 


This  result  shows  that  for  the  SH  wave  equation,  and 
consequently  for  liquid  media,  the  error  truncation  order 
is  independent  of  the  boundaries  (in  the  case  of  stress 
and  displacement  continuity).  Then,  more  generally, 

N 

e  =  l  0  (h  n,At  )  (3.9.12) 

n=  1 

From  a  Taylor  series  (3.9.12)  can  be  written 

£  =  f  [D2nVAt2(n'1)  -  H.  [D2nV  +  D2nv]h2(n-1)] 

L  _  1  t  p  x  Z 

n=  1 


, .  2n  2n 
+  0  ( At  , h  ) 


(3.9.13) 
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But  from  the  equation  of  motion 


DJnv  =  [D2V  +  D2V](n)(-)n 
t  x  z  J  n 


then 


N 


H.  /iLv  n“  1  ^  v,2,,!  (n)  <n-D 


£=  IJ 


-  [D2nV  +  D2nV]h2(n-1)] 

X  z 


(3.9.14) 


which  lead  for  an  undirect ional  wave 


N 


E  = 


n=  1 


I  t  h 2 n _  1  [  ( — )  n -  1  2 (n"1)  .  1]D2nv 

P  p  h  - 


(3.9.15) 


Then,  in  the  case  of  an  unidirectional  wave  e  =  0 
for  any  order  of  development,  if 


77  <^>2  =  1 
P  h 


(3.9.16) 


This  relation  is  important  since  it  forecasts  the 
identity  between  the  analytical  solution  and  the  numerical 
solution  for  synthetic  seismograms.  This  fact  can  be 


easily  verified. 
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Then  in  the  case  of 


y 

p 


<^>2 


1 


we  have  not  only  consistency,  we  have  an  identity  between 
unidirectional  SH  waves  and  their  analytical  solution. 
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3.10  CONSISTENCY,  STABILITY  AND  CONVERGENCE 

When  approximating  the  analytical  problem'  to  a  discrete 
problem  one  main  requirement  is  when  h^  ,  At  -*  0/  AU  -►  AU 
i.e.  the  discrete  problem  becomes  equivalent  to  the  con¬ 
tinuous  one,  and  the  numerical  problem  is  said  to  be 
consistent  to  the  analytic  problem. 

However,  the  consistency  does  not  necessarily 
imply  that  the  numerical  solution  approximate  the  analyti¬ 
cal  solution  and  converged  to  it. 

Although  theoretically,  consistency  and  convergence 
are  sufficient  conditions.  Van  Howen  (1968)^;  rounding 
errors  give  rise  to  a  computational  solution  instead  of 
the  time  difference  solution.  This  leads  to  the  stability 
of  a  difference  scheme. 

3. 10. a  Consistency 

The  consistency  is  measured  by  the  truncation  error. 

% 

Since  F.D(U)  =  0  (by  definition  (3.9.6)) 

E  =  |f.D(U)  -  FD  ( U )  |  (3.10.1) 

The  equations  (3.9.7),  (3.9.8)  shows  that  if  h, 

At  ->  0  £  ■+  0 
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then 


F  D (U)  +  FD (U) 


1 .  e  . 


F  D  ( U )  -*  0 


and 


F  D  (U)  -►  (A  +  D  )  U  =  0 


(3.10.2) 


Hence,  the  finite  difference  development  of  the 
equations  of  motion  in  elastic  media  trend  to  the 
differential  equation  when  h,  At  0. 

3.10. b  Stability 

There  are  two  methods  of  investigating  the  stability 
criteria.  One  uses  a  finite  Fourier  series  (method  of 
Von  Neumann); the  other  expresses  the  equations  in  matrix 
form  and  examines  the  eigenvalues  of  an  associated  matrix. 
We  shall  use  the  last  method  since  it  takes  into  account 
the  boundary  conditions  but  first  we  shall  examine  the 
simplest  problem. 

3.10. b.l  Stability  of  two  time  level  equations 


We  have 


l+l 


l 


U 


AU 


(3.10.3) 


• 

. 

' 
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0 

U  being  the  initial  condition 


U1  =  AU° 


U*+1  =  fl£+1u° 


The  system  is  stable  if 


(3.10.4) 


1 1*1 1  <  i 


3 . 10 . b  .  2 

The  system  can  be  written  on  the  form: 


U 


l+l 


^  l 

AU  + 


% 

BU 


l-l 


(3.10.5) 


% 

(since  A  =  F  D(A)). 

By  analogy  with  (3.10.3)  we  consider  a  vector  V 
such  that 


U 


U 


l-l 


(3.10.6) 


Then  (3.10.3)  can  be  written 
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U 


l+l 


(3.10.7) 


where  I  is  the  unitary  matrix  and  the  amplification 
matrix  is 


(C)  = 


A  B  \ 


l  I  0 


(3.10.8) 


By  analogy  with  (3.10.4),  (3.10.7)  yield 


U 


1+ 1 


I  A  B 


I  0 


l+l 


V 


0 


(3.10.9) 


T-  0  0  .  . 

where  V  =  U  (initial  conditions) . 

The  boundedness  of  (3.10.9)  implies  (see  Boore  1970) 


C  <  1 


l .  e  . 


A  <  1 


(3.10.10) 


where  A  represents  the  eigenvalues  of  (C) ,  i.e. 


A  is  solution  of 


det | (C)  -  Al |  =  0 


(3.10.11) 


From  Smith  (1978)  we  introduce  the  matrix 


V 
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(3.10.12) 


Since  det  (D)  /  0,  we  write 


det | ( (C)  -  Xl)  •  (D)  =  0 


(3.10.13) 


which  lead  to 


det | A I  -  (A)  -  (B)/a|  =  0 


(3.10.14) 


From  the  explicit  form  of  the  equation  of  motion  we  have 


(B)  =  -I 


then  (3.10.14)  yield 


det  (A) 


0 


(3.10.15) 


with 


A  +  1/A 


then  if  A  is  an  eigenvalue  of  (C)  so  is  1/A.  Consequently 


1 


(3.10.16) 
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and  the  stability  condition  (3.10.10)  can  be  written 


A  <  2 


(3.10.17) 


The  eigenvalues  of  A  are  given  by  the  Gershgorin  and 

Brauer  theorem  which  states  that  if  the  radius  p.  is 

3 


as  follows : 


n 


p3  =  J,  'v|(1  -  V 


(3.10.18) 


then  the  eigenvalues  of  (A)  lie  in  at  least  one  of  the  circles 


z  -  a  .  .  <  p  . 

3  3  1  ~  3 


(3.10.19) 


% 

(A  +  2I)U  where  A  is  the  F.D.  of  the 


Then 


,  2  20j  2r\j 

'At  (a  D  +  8  d  )  +  2 

XX  zz 


2% 

(A  +  y )  At  D 


xz 


AU  = 


U 


i  (  A  +  y  )  A  t  D 


xz 


9  9  rv  9  r^' 

2+At(aD  +  8  D  )/ 

ZZ  XX, 


(3.10.20) 


where 


' 
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a 


2 


(X  +  \1)/1 


e2  =  v/i 


With  the  Gershgorin  theorem  we  obtain 


2 


(a2  +  e2)|  < 


<  2  . 


(3.10.21) 


Then  the  stability  condition  is  given  by 


1 


max ( a  + 


e2) 


(3.10.22) 


or 


aAt 

-  < 

h  - 


a 


(3.10.23) 


We  note  that  in  the  case  of  SH  wave  the  stability  condition 
becomes 


At 

h 


< 


(3.10.24) 


where  n  is  the  number  of  spatial  dimensions  n 


1,  2,  3 
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3.10.C  Convergence 

There  is  convergence  if  the  solution  of  the  finite 
development  tends  to  the  solution  of  the  differential 
equation  when  At,  h  ■+  0 . 
i  .  e  . 


6  =  U  -  U  I  0  when  At  ,  h  ->  0 

e  i  i 

For  the  wave  equation  case  the  discretization  error  is 
given  by 


6  =  At 2e 

e 


where  e  is  the  truncation  error. 

We  note  that  the  finite  development  of  the 
equation  is  convergent  since  stability  and  consistency 


implies  convergence  (Lax  Equivalence  Theorem;  Richtmeyer 


and  Morton  (1967) ) . 


. 
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3.11  ALIASING  AND  GRID  DISPERSION 

3. 11. a  Discretization  of  the  spatial  waveform  implies 
an  "effective"  sampling  at  the  node.  The  condition  is 
therefore  that  the  frequency  domain  where  the  amplitude 
spectrum  of  a  function  F  is  non  zero  is  finite  (i.e.  F 
has  to  be  band  limited) . 

The  relationship  between  the  temporal  or  spatial  grid  and 
the  frequency  or  wave  number  is  given  by  the  sampling 
theorem  (Bracewell,  1965;  Kanasewich,  1975)  which  states 
that 


where 


FW 


K.} 

1 


0  if  I K . I  >  K 
1  l  1  c 


(3.11.1) 


It  is  therefore  recommended  when  we  have  a  source  with  a 


broad  spectrum,  to  reduce  the  effective  frequency  domain  by 

sin  x  sin  x  1 

convolving  by  a  function  of  the  form  -  or  ( - )  , 

X  X 

the  amplitude  spectra  of  these  functions  being  the  well 

known  box  car  or  triangle  (Bracewell,  1965) . 

Since,  in  seismic  prospecting  the  incidence  angle  is 
close  to  the  vertical,  the  apparent  wave  number  along  the 
x  axis  is  equal  to 


K 


1 


K  sin  6 


. 

' 


9  5 


which  implies  a  wide  grid  spacing  in  the  x  direction  if 
the  economic  consideration  (computer  time)  surpasses 
the  interest  in  surface  waves. 

Generally,  as  a  rule  of  thumb,  eight  to  ten  samples 
per  wave  length  are  considered  minimum  to  approximate 
the  wave  field  (Boore,  1972;  Claerbout,  1976) . 

3.11.b  Dispersion 

The  dispersion  appears  as  a  variation  of  phase 

velocity  and  group  velocity  as  function  of  the  frequency. 

Consequently,  while  the  low  frequency  part  of  the 

signal  will  travel  with  the  velocity  v  =  V  =  V,  the 
3  0  g 

under  sampling  will  affect  the  high  frequency  part  of  a 
signal  whose  group  and  phase  velocity  will  be  V  <  V  < 
VQ .  Then  an  apparent  attenuation  will  be  introduced 
as  well  as  a  delay  in  phase  resulting  in  a  tailing 
effect. 

To  investigate  these  effects  we  will  follow 
Brillouin  (1953)  and  Alford  et  al  (1974)  in  the  case 
for  SH  waves. 

3 . 1 1  .  b . 1  SH  wave  dispersion 


Let  us 


consider  a  harmonic  plane  wave  of  the  form 


■ 
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U  =  Uq£1  (.Wt-kjX-k^  z ) 


(3.11.2) 


where 


=  K  cos  6 


K  ^  =  K  sin  6 


(3.11.3) 


The  propagation  is  governed  by 


u  _2 
— T  =  V  U 


V 


0 


(3.11.4) 


Second  order  finite  difference  development  of  the  above 
equations  lead  to 


l-l  „  l+l  U  -  2U  +  U 

U _ -  2U  +  U  _  m-1 _ m+1 

A  2  "  2 
V0At  h 


U  ,  -  2U  +  U 

n-1  n+1 


(3.11.5) 


where  U(x,z,t)  =  U(mn,  nh ,  ZAt) 

To  simplify  the  notation  U(m,n,£)  is  represented  by 

U  and  U.  ..  « .  =  U  .  ,  etc. 

(m- 1 , n , t)  m-1 


If  P  = 


V  At 
0 


we 


obtain,  with  P  <  /2/2  (stability  con¬ 


dition)  ,  after  substituting  (3.11.2)  into  (3.11.5) 
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sin 


2  uJ  At 
2 


K  h 
1 

2 


+ 


.  2 
s  m 


K2h 


(3.11.6) 


If  G  is  the  number  of  samples  per  wave  length  V  and  V 
the  phase  and  group  velocities,  we  can  write 


9 


Kh/2  =  P/G  • 


Then 


.  uJ 

v/vo  ■  M 


.  "1  r„  r  •  “2  (Ttcos  6)  ,  .  2  (fT'sin  0  )  -,  •, 

sm  {p[sin  - — -  +  sm  - ]} 


0 


...(3.11.7) 


and  the  group  velocity  will  be  obtained  by  differentiating 

(3.11.6)  with  respect  to  K. 

V  _  _ 

— =  [sin  (-  cos  0)  cos  (r-  cos  0)  cos  0 

vo  G  G 

+  sin(|A  sin  0)  cos  (A  sin  0)sin  0]/ 

G  G 

r  r,  2  .  2  /n  Q  ,  2  .  2  .  rr  .  2  ,  h 

{ [ 1  —  P  sin  (—  cos  0)  -  P  sin  (—  sin  0) J. 

G 

[ sin2 (—  cos  0)  +  sin2  sin  0  )  ]  **  }  (3.11.8) 

G  G 

If  the  propagation  is  parallel  to  the  grid  0=0  and 

(3.11.7) ,  (3.11.8)  yield : 


■ 
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VV0  =  f“  si"_1[P  sin  2] 


(3.11.9) 


V  /V  = 
G  0 


cos  — 
G 


(1  -  P 


.  -  2  77  h 

sln  G 


(3.11.10) 


3.7.b.2  P . SV  Dispersion 

In  the  absence  of  body  forces  (2.1.1)  yield: 


'v.  -> 

AU+pD^U  =  0 


(  3 . 1 1 J  1 1 ) 


with 


A 


/  ,  ^  ^ 

/ ( X  +  2y ) D  +uD 

XX  zz 


( A+y ) D 


xz 


( A  +  y ) D 

xz 


\ 


j 

%  ^  / 

(  A  +  2  y  )  D  +  y  D 

Z  Z  XX/ 


(3.11.12) 


r\j 

D  ,  D  ,  D  being  the  finite  difference  operator  of  the 
xx  zz  xz 

corresponding  derivatives. 

As  in  section  2.7.b.l  we  can  write  for  a  harmonic 
plane  wave  with  Ax  =  Az  =  h  =  constant 


- 
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'Xi 

D 

xz 


r\j 

D 

xz 


lx”  =  ^2  [Ux  -  2U  +  03] 


D 


xx 


=  2 [ 1  -  COS  K  h ] /h' 


4.2  h 

^Sin  K1  I 

h 


D 


z  z 


4  •  2  v  h 

"2  Sin  K2  2 
h 


^  „  .  2  wA  t  .  .  2 

Dtt  =  4  sin  — /At 


2h 


—  [cos [ K x h  +  K2h]  -  cos[Klh  -  K2h]] 


sin 

— —  (K^h) sin (K2h) 


for  convenience  we  write  (3.11.12) 


% 

A 


% 

a 


(3.11.13) 


(3.11.14) 


(3.11.15) 


(3.11.16) 


(3.11.17) 


with 


(3.11.18) 
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(3.11.11)  has  the  form 


r\r+  ->■ 

AU  -  AU  =  0  (3.11.19) 


or 


/A‘a  ii  ai2\ 

\a21  A-a22/ 


(3.  11. 19)  1 


which 


has 


equation 


2  n,  %  %  %  %2 

A  '  A(all  +  a22>  +  aHa22  -  al2  "  ° 


(3.11.20) 


The  roots  being  given  by 


A 


(ail  +  a  2  2  +  A 


(3.11.21) 


A 


<\>  <\j  r- 

(ail  +  a  2  2  >  '  A 


(3.11.22) 


where  the  discriminant 


■  (X  +  y)2[(Sxx  - 


a-2  A/ 

D  )  +  4  ( D 

zz  xz 


2 


] 


A 


(3.11.23) 
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if  we  write 


S1  =  sin  K^h/2 


S ^  =  sin  h/2 


(3.11.23)  yield 


4  =  iA±|i_[(4s2  .  4S2)2  + 


4  (  s  i  S  2 ’  2 [1  +  S2S2 


(S2  +  S2,] 


or 


A  =  42  — ±M.j—  [  ( g  ^  .  s2,2  +  4(S1S2)2[1  +  S2S2  -  (S2  +  S2)] 


then 


2  (  X  +  U )  r  2  2  2  2  %  ,  2 

A  =  4  '  y  [ (Sx  +  S2  ±  2S1S2)] 


and 


^  -  ±4  ^  S2  *  2S2S2]2 


(3.11. 24) 


then 


2  _  2  .  _  .  r  _  2  _  2  _ 2  _  2  n  ,  2 


A  =  { 2  ( X  +  3y)  (S1  +  S2)  +  2  ( X  +  y)  [si  +  S2  -  2S1S2]}/h' 


' 
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A1  -  (4(X  +  2y )  ( S  2  +  S2)  -  4  (  X  +  y)S1S2)/h2  (  3  .  1 1 . 2  5  )  a 

2 

a2  =  ( 4y  ( S  x  +  S  2 )  +  4 ( X  +  y)S1S2)/h2  (3.11.25)b 


<\i  <\j 

But  pD^^U,  pDttW  are  eigenvalues  and 


^42 
PD  U  =  - -  p  sin  (up  At/2) 

At  1 


pD  W  =  sin 2  (u>  At/2)  (3.11.26) 

At  2 


then  (3.11.25)  and  (3.11.26)  yield 

sin2  (u^  At/2)  =  ( X  +  2y ) At 2/ph2 [ S 2  +  S2]  -  [ X  +  y ] At 2 /ph 2 S 

...  (3.11.27) 

sin2  (up  2  At/2)  =  yAt2/ph2  +  (A  +  y)At2/ph2  S  ^  S  2 .  (3.11.28) 


Recall  that 


Kxh/2 

K2h/2 


e 


—sin  6  , 

G 


(3.11.29) 


' 
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As  usual  ,  we  set 


a 


2 


(X  +  2y)/p 


3 2  =  y/p 


f 


Then  the  ratio  between  the  numerical  phase  velocity  and 
the  media  velocities  are: 


ap/a 


0 


(3.11.30) 


6P/f5o 


U) 


2 

k|$ 


0 


(3.11.31) 


This  choice  is  a  consequence  of  the  fact  that 


A 


-  =  4  ( X  +  y)  (S 


-  S2> 


i  .  e  . 


X 


2 


Taking  into  account  (3.11.29)  the  expressions  (3.11.30), 
(3.11.31)  can  be  written: 


ap/aQ  =  oj  1hG/2TTao 


(3.11.32) 


6p/e0  =  a)2hG/2rrBo 


(3.11. 33) 
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•  .  rU 

Since,  as  mentioned  above  D  U  and  D  W  are  eigenvalue 

tt  tt  L 

a  At 

equation  (3.11.25)  yields,  after  setting  P  =  0 


h  ' 


Y  =  3/a 


sin2(d)1  At/2)  =  p2[s2  +  S2]  -  [l  -  y2]p2S1S 


(3.11.34) 


sin2(U)2  At/2)  =  y  2p  2  [  s 2  +  S2]  +  [l  -  y2]?^^.  (3.11.35) 


T  hen 


to  .  = 


sin  1p[(s2  +  s2)  -  (l  -  y2)s1s2]35 


(3.11.36) 


= 


2  .  -1  r  /  2  .2,  2  , ,  2  -i  ^5 

—  sin  p[(s1  +  s 2 ) y  +  (l  -  y  )s1s2] 


(3.11.37) 


The  dimensionless  phase  velocity  ratio  is  then  given  by 


ap/aQ  =  sin  1p[(Sp  +  S22)  -  (1  -  Y2)S1S2]1  (3.11.38) 


V6o 


r  -  1  2  2  2 

’  sin  ■Lp[y^(S1  +  S  2 )  +  (1 


pit  y 


Y2)siS2]i5 


(3.11.39) 


and  the  group  velocity  ratio  is  given  by 


a  /a _  = 


i  3ui 


'G  0  aQ  3k 


(3.11.40) 


o\$  * 


■ 
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VB0  =  r0JT-  (3.11.41) 

We  note  that  for  y  =  1  or  0  =  0  the  results  are  the  same 
as  for  the  SH  wave  or  the  scalar  wave  equation  given  by 
Alford  et  al ,  the  dispersion  being  minimal  for  values  of 
P  chosen  at  the  stability  limit. 


■ 
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0  =  0 


Fig.  3.5  Dispersion  curves 

3. 5. a  Normalised  phase  and  group  velocity  for 
a  propagation  parallel  to  the  grid  0=0 
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P  =  ihAL  -  „  7 

h 

3.5.b  Normalised  phase  and  group  velocity 
different  angle  of  incidence 


at 


CHAPTER  IV 


Comparison  with  Recursive  Developments  Derived 
from  Variational  Methods  -  Case  of 
Potential  Equations  of  Motion 

We  have  mentioned  in  Chapter  II,  three  main  methods 
of  numerical  analysis: 

-  The  direct  expansion  by  a  Taylor  series 

-  The  Rayleigh  -Ritz  method  (application  of  Dirichlet  principle)  . 

-  The  Galerkin  method  (usually  called  the  variational  method 
After  showing  the  equivalence  of  the  two  last  methods 

we  shall  develop  in  this  chapter  a  finite  development  derived 
from  variational  principles.  We  will  then  compare  it 
with  the  finite  difference  scheme  of  Chapter  II  obtained 
by  direct  application  of  Taylor  expansion  on  the  usual 
formulation . 

Friedrich  (1962)  ,  Friedrich  and  Keller  (1966) 

2 

considered  the  problem  for  the  operator, -V  ,  and  proceeded 
by  numerical  integration  in  triangular  elements. 

The  approach  used  in  these  pages  will  be  more 
operational,  in  such  a  way  that  we  take  full  advantage 
of  the  bilinear  forms  obtained  by  use  of  the  z  transform 
and  its  conjugate. 

Use  of  Sobolev  space,  W^  will  give  great  advantage 
for  the  device  of  a  trial  function. 
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■ 

’ 


. 
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In  the  second  part  of  this  chapter,  we  will  consider 
the  equations  of  motion  (potential)  and  compare  their 
finite  development.  The  case  of  the  boundary,  considered 
as  a  limit  of  a  transition  zone,  will  also  be  considered 
for  those  equations. 

As  in  Chapter  II,  we  will  consider  the  equation  of 
motion  in  the  form 


AU  ( P )  =  F  ( P ) 


(4.1.1) 


where  A  is  positive  definite 

U  £  W ^  (Sobolev  space) 

2 

P  £  ft  c  IR  x  [o  ,  T) 

According  to  the  Ritz  method,  we  have  to  minimize 
the  functional 

I  (  U )  =  ( AU  ,  U )  -  2  (  U  ,  F  ) 

=  /  [ U ( P ) AU  -  2U  (P) F  (P)  ] dft  (4.1.2) 

n 

U ( p )  represents  the  acceptable  functions  from  the 
field  of  definition  D  of  the  operator  A. 


Since  for 


' 
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I(U)  =  min [ ( AU , U )  -  2(U,F)], 


(4.1.3) 


U  is  the  solution  of  4.1.1. 

Then,  the  problem  (4.1.2)  can  be  written: 


(  AU , U )  -  (  F  ,  U )  =  0 


(4.1.4) 


which  is  the  Galerkin  form. 

4.1  FINITE  DEVELOPMENT  OF  THE  EQUATIONS  OF  MOTION 
DISPLACEMENT  BY  VARIATIONAL  PRINCIPLE 

4.1. a  Dirich  let  Integral  Form 

After  substituting  the  elasticity  operator  by  its 
value,  (4.1.4)  yields: 

-/  {d  [ [ A+2y ] d  u  +  Ad  w]u  +  d  [y d  u  +  yD  w]u}dft 

X  X  Z  Z  Z  X 


-/  F . Udft  =0  (4.1.5) 

ft 


{d  [  [A  +  2 y ] d  w  +  Ad  u] u  +  d  [yD  w  +  yD  u]u}dft 

Z  Z  X  XX  z 


-/  F.Wdft  =  0 

ft 


(4.1.6) 


' 


■ 


Ill 


By  Gauss  Theorem  (4.1.5)  and  (4.1.6)  yields 


/  { ( A+2 y ) (d  u) 2 
ft  x 

o 

+  Ad  WD  U  +  y(D  U)  +  yD  WD  U>Udft 

Z  X  Z  X  z 

-  /  [(A+2y)D  u  +  Ad  wlu.dT 

Y  x  z 

-  /  y  [d  U  +  D  W ] U . d  T  -/  F.Udft  =  0  (4.1.7) 

r  z  x  ft 

/  {  (  A  +  2y)  (D  w) 2 
ft  x 

+  Ad  ud  w  +  y(D  w) 2  +  yD  ud  w}wdft 

X  Z  X  Z  X 

-  /  [(A+2y)D  w  +  Ad  ulw.dT 

p  Z  X 

-  /  y[D  W  +  D  UjW.DT  -  /  F.Wdft  =  0  (4.1.8) 

r  x  z  ft 

If  we  express  the  continuity  of  displacement  and 
stresses  at  the  boundary  (Boundary  conditions)  the 


integrants  of  / 

r 

(  .  ) dT  become  s 

/  (  . ) dT  =  0  (4.1.9) 

r 

and  (4.1.7)  and  (4.1.8)  yields: 
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/  {(X  +  2y)  (d  u)2  +  Ad  wd  u  +  y(D  u) 

X  Z  X  Z 


Z  X  z 


+  yD  WD  uludft 


X  z 


/  F.Udft  =  0 

Q 


(4.1.10) 


/  {(A+2y) (d  w)  +  Ad  UD  w  +  y(D  w) 

n  Z  X  Z  X 


+  yD  WD  u}udf2 


z  x 


/  F.Udft  =  0 

Q 


(4.1.11) 


4.1.b  Discretization  and  Choice  of  a  Trial  Function 

Since  the  objective  is  the  comparison  with  the 
finite  difference  development,  it  is  normal  that  we 
choose  a  trial  function  from  the  same  field  of  definition  . 
Then  we  have 


U  £  W ^  (Sobolev  space) 


(4.1.12) 


which  implies  the  conservation  of  norm 


(4.1.13) 


Then  (4.1.12)  implies  that  the  function  and  its  first 
derivative  has  to  have  respectively  the  same  expression 


as 


the  finite  difference  expression. 
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As  in  Chapter  III,  we  consider  an  elementary 
region 


such  that 


ID  = 


Then,  for  instance. 


D  U 
x 


] 


for  U  e  ID 


11  ii  111  ,  „  , 

Up  e  t  c  .  (  4  .  ] 


If  we  consider  U  such  that 


i  (  wt  -  K  rnhj  -  K  nh  ) 

u  -  V 


we  can  use  the  z  transform 


with 


z 


e 


-  i  K  h 
2 


2 


x  = 


‘iKihi 


iWAt 

T  =  e 


.14) 
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then 


U 


m  ,  n-1 


=  UZ 


U 


=  UX 


m-1  ,  n 

£+1 

U  =  UT 
m  ,  n 


Then  (4.1.14)  can  be  written 


d  u  =  ^ —  ( 1  -  z)  for  u  e  id  1 1  Uid111 

X  h 


(4.1.15) 


and 


UZ 

D  U  =  - 

x  h 


-1 


(l  -  Z)  for  U  e  Dx  U  D 1 V  etc. 


(4.1.16) 


For  more  simplicity  we  will  have  h  =  h  ,  h  =  h  , 

13  2  4 

without  loosing  the  generality  of  our  problem. 


2  2 

4.1.c  Bilinear  form  of  the  expression  (D  U)  ,  (DU)  , 

x  z 

D  WD  U,  D  WD  U 

X  Z  Z  X 


Using  the  z  transform,  we  can  write 


(D  U) 
x 


4- 

hi 


-  X][l  -  X] 


U  E  |D 


(4.1.17) 


then 


. 
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(DxU)2 


[(1  -  x]uQ 

1  I  I  ,  III 
e©  U© 


X*[l  -  X]UQ] 


EIDI'iU\D 


I 


(4.1.18) 


and  according  to  section  (3.2)  if  we  consider  the 
parameter  p  for  instance 


V ( D  u)2  =  [p  [l  -  X ] U  -  p  X*  [l  -  xlu  ]  (4.1.19) 

x  ,2  3  0  1  0 

hi 

when  y  =  (V1  +  P  ) / 2 

u3  =  <y;  +  vp/2 

then 

P  (D^U)  =  "  \  "  U0^  "  P3^U3  ”  U0^  (4.1.20) 

h. 


By  the  same  way 


and 


U 


(D  U) 
z 


”  [i  -  z]  [i  -  z] 


(4.1.21) 


P  (D  U) 
z 


U 


J  [  (1  -  Z'lVo  -  Z*t1  •  Zly4U0] 


(4.1.22) 
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then 


U 


y(D  u) 
z 


2  [u2[u2  -  V  -  y4[U4  -  °0” 


To  calculate  the  expression  D  UD  U  we  have  to  note 

x  z 


ID  (D  UD  U)  1  =  ( X  ,  Z )  PI  D 1 

x  z 


which  implies 


ID  (D  U)1  c  (x  f  Z)  (1  D1 
x 


and 


|D(D  U)  c  (x,Z)  fl  D' 


then 


U 

D  U  =  [ ( 1  ~  X )  +  Z[1  -  X] ] 

X  z 


then 


U 

D  U  =  ^  [[1  -  Z)  +  X(1  -  Z]] 

z  2 


U  W 
0  0 

D  WD  U  =  — r — 7—  [  [ 1  -  x]  +  2(1  -  x]][(i  -  z)  +  x ( 1 
x  z  4hi  2 


(4.1.23) 


that 


(4.1.24) 


(4.1.25) 


-  Z)  ]  * 


...(4.1.26) 
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DxWDzu  =  4h°h  1(1  '  x]  +  Z<1  -  x3 3  r ( 1  +  X 


and  referring  to  section  3.2 


W  U 

yD  WD  U  =  ■  [(1  -  X]  +  Z ( 1  -  X] [ [1  +  x' 

x  z  4h ri  L 

1  2 


By  expliciting  (4.1.28)  we  obtain 


U 


yDxWDzU  =  {[[W6  -  w3]  +  [W1 


+  [ (W2  -  w?)  +  (WQ  -  w3) ]u3 


-  [<wl  -  wQ)  +  <w5  -  w4>]y; 


-  1(W0  -  w3)  +  (W4  -  w8)]v4> 


In  the  same  way  we  have 


UW 

D  WD  U  =  - —  [d  -  Z)  +  X(1  -  Z ) ] L ( 1  ~  X] 

z  x  4h  yh  2 


]  -  z*[i  +  x  1]] 

...(4.1.27) 


1 J y 4  -  z*[i  +  x  1) ] 

...(4.1.28) 


w0]]y; 


(4.1.29) 


+  Z  [  1  -  X]  ]  * 

(4.1.30) 


and 
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VDzWDxU  “  4Th  l(1  -  2>  +  X[1  -  Z][[l  +  Z-1], 

12  3 

-  X*  [l  +  z  1 ] y x ]  (4.1. 31) 

Explicitly  (4.1.31)  yields 

yD  wd  u  =  -  — - -  ■■  {[(W  -  w  )  +  [w  -  u  lly  + 

z  x  4hlh2  6  1  2  o J  J  M1 

+  [ (W2  -  wo)  +  (W7  -  w3)]y^ 

-  [ (W1  -  w5)  +  ( w 0  -  w4 ) ] y ~ 

-  [<W0  -  W4)  +  (W3  -  w8 )  ]  y  3 }.  (4.1.32) 


As  an  abbreviation  we  will  call  U.F V  (  •  )  the 
development  obtained  by  finite  variational  method,  in 
this  chapter. 

We  remark  from  equation  (4.1.28)  and  (4.1.32)  that 


FV [ y D  WD  U]  =  F.V(yD  WD  u] 
Z  X  x  z 


vp  e  n 


but  p  t  r 
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Comparing  with  the  results  obtained  in  Chapter  III,  we 
have 


and 


if 


uf.d{£d  u}  =  -ufv{£d  u} 

xx  xx 


UF  D{£D  u}  =  -UFV{£D  U} 
z  z  ZZ 


Vp  e  ft  =  (0,  +  r ) 


UFD{£D  w}  =  -UFV { WD  u} 
xz  X  z 


Vp  e  ft 


p  t 


ufd{ £d 

xz 


p  e 


r . ( n . Oz ) 


W}  ^  -UFV { WD  U} 

X  z 

r . ( n . Oz ) 


(4.1.33) 


(4.1.34) 


(4.1.35) 


when  n.Oz  design  the  projection  of  the  normal  of  the 
boundary  to  the  Oz  axis. 

In  this  case  if  hi  Oz  then  P  can  be  considered  as 
an  inner  point. 

Equations  (4.1.33)  ,  (4.1.34)  ,  (4.1.35)  can  ‘be 


Written : 


•  ■  1 
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U  FD{£D  U}  =  -UF  V { £ ( D  U)2} 

XX  X 

U  FD{£D  U}  =  -UF  V  {  £  (  D  U)2} 

Z  2  2 

VP  e  ^  =  (fi  +  T)  (4.1. 36) 

U  fd{£d  W}  =  -UF  vUd  wd  u} 

X2  X  Z 

VP  £  (P  inner  point) 

and  P  t  T  (4.1.37) 

u  fp{£d  w}  l  -uf.v{£d  WD  u} 

XZ  X  z 

VP  £  r  (P  Boundary  point)  (4.1.38) 

2 

since  we  have  c  ir  x  [0T] 

in  the  absence  of  body  forces  on( can  write: 

-  /  F  XS  dft  =  /  pD  U  U  dft  (4.1.39) 


i  .  e  . 
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-  /  Fdft 


/  P  (DU)  dft, 
O  r 


TFu0tul 


-  2U  + 


u"1] 


{Mass } 


C2, 


(4.1 


4.1.d  Recursive  Developments 

Replacing  in  equations  (4.1.10),  (4.1.11)  the 

respective  expressions  by  their  FV{.)  developments  and 
taking  into  account  (4.1.10) ,  we  obtain 

A  2 

U1  =  2U  -  U_1  +  -y  -  - - [  [U1  -  UQ  ]  ( X  +  2 y )  x 

hl E  Pl  +  P3 1 

'  (U0  ”  U  3 }  (X  +  2^i)  31 


2  At' 


h2[P2+P4l 


[<U2  -  u0]u2  -  [u0  -  u4)]y4i 


- -^7 - -  f'vUd  wd  u} 

4h1h2(p3+pi)  z  x 


At 


2 


f  . y { y D  wd  u} 

x  z 


.40) 


+ 


4hlh2(P2+P4) 


(4.1.41) 
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+ 


2  At' 


h  9  (P  +P J 


[  (W. 


WQ)  (A  +  2y)2 


-  (wo  -  V  <*+2y>4] 


2  At 


hl (P1+P3) 


[(wi  -  W0]V1  *  <wo 


-  w3)P3] 


+ 


4hlh2  [P4  +  P21 


F  v[Ad  ud  w] 

X  z 


+ 


At 


4hlh2[p3  +  PJ 


F  V [ p  D  UD  w] 
Z  X  J 


(4.1.42) 


By  comparison  with  equations  (3.4.1)  and  (3.4.2) 
we  see  that  the  two  systems  are  identical. 

4.1.e  Uniqueness 

As  a  consequence  of  the  identity  of  the  development 

between  the  variational  method  and  by  Taylor  expressions  we 

are  led  to  the  following  proposition: 

The  elastic  wave  equation  admits  an  identical  second 

order  numerical  expression  if  the  trial  function  has 

* 

the  same  norm  in  the  conservation  of  energy  space  W^. 
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4.2  CASE  OF  EQUATIONS  OF  COMP RE S S I ONAL  AND  ROTATIONAL 
DEFORMATIONS 

Let  us  write  the  equations  of  motion  for  a  homogeneous, 
isotropic  elastic  media 

(X  +  2y )  VVU  -  pVAVAU  =  pU  -  pK  (4.2.1) 

If  we  take  successively  the  divergence  and  the  curl  of 
equation  (4.2.1) ,  the  dilatation 

P  =  div  U  (  4 . 2 . 2 ) b 


and  the  vector  rotation 


S  =  curl  U  (  4 . 2  .  2 ) c 


satisfy  the  respective  equations 

(X  +  2y)V2P  =  pP  -  pV.K  (4.2.2) 

yV2S  =  pS  -  p V  A  K  (4.2.3) 


where  K  represents  the  body  forces. 


■ 
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The  aim  of  this  section  is  to  give  a  finite  develop¬ 
ment  for  those  equations  by  the  two  methods  developed  in 
these  pages.  i.e.  The  finite  development  by  variational 
principle,  and  the  generalised  solution  of  the  direct 
Taylor  development. 

4. 2. a  Development  by  Variational  Method 

If  we  consider  the  usual  elementary  domain  D  =  , 

the  respective  scalar  product  of  (4.2.2) ,  (4.2.3)  with  P 

and  S  yields  in  the  absence  of  body  forces 


/  a2  ( VP)  2d£2 
o  h 

^ h 


a2Vpdf 


(p)  d n 


(4.2.4) 


/  e2(vs)2dft 


32VsdT 


/  (s)2dft 

Q  n 


s , p~ e  wj  c  l2 (ft) 


with 


where 


ft,  =  ft,.  +  r.  c  ft  c  R  x  [  OT  ] 

h  h  h 


3  =  V/P 


(4.2.5) 


2 

a 


(A  +  2]a)  /  2 


' 

. 
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(4.2.4),  (4.2.5)  yield  from  section  4.1 


/ 

ft. 


pf  via2 ( Vp) 2 } 


-  /  Pa  VPdT  =  /  PF  V  [d  p]2dfi  (4.2.6) 

r,  h  a  t  h 

n  h 


/  sf  v{32 (Vs) 2} 

ft. 


-  /  se2vsdr,  =  / 


fa 


SF  v[d  s]  da 

t  h 


in  ft  c  R  x  [OT] 


(4.2.7) 


The  boundary  conditions  of  (4.2.6)  can  be  written 


/  a2VPdT 

r 


Ath2[“lPx  '  “3Px]  +  hlta2Pz 


a2P  ] At 
4  z 


(4.2.8) 


where  a  ,  a^  represents  the  values  of  the  parameters  in  the 
usual  domains,  i.e.  D1  U  DIV,  etc. 

Since 


P  =  U  +  W 
x  z 


(4.2.9) 


(4.2.10) 


The  stress  continuity  conditions  can  be  written 
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( ot 2 p  -  262w  ).  =  (a2p  -  232w  ) 

z  1  z  3 

(a2p  -  232u  )  0  =  (ot2p  -  2B2u  )  „ 

x  2  x  4 

[6  (S  -  2UZ) ]2  =  [B(S  -  2  u  z ) ]4 

[ 3 ( s  -  2W  ) ]  =  [3(S  -  2W  ) ] 

X  X  j 


(4.2.11) 


(4.2.12) 


(4.2.13) 


(4.2.14) 


The  expressions  (4.2.13)  and  (4.2.14)  being  identical. 

After  differentiating  (4.2.11)  and  (4.2.14)  respectively 
with  regard  to  the  variable  x  and  z ,  we  obtain 

a2p  -  a2p  =  B2s  -  32s  .  (4.2.15) 

lx  3x  3z  lz 


In  the  same  way  taking  into  account  (4.2.11)  and  (4.2.14) 


2  2  2  2 

c^p  -  a,p  =  3.  s  -  32s 

2z  4z  2.x  4x 


(4.2.16) 


the  boundaries  concerning  equation  (4.2.7)  yields 


x 


$2s 

3  x 


(4.2.17) 


S 

z 


34S 
4  z 


(4.2.18) 


. 


127 


with 


26 


2 


After  inserting  the  boundary  conditions  in  the  equations 
(4.2.6),  (4.2.7)  and  replacing  the  values  of FV( • }  by  their 

developments  (section  4.1)  we  obtain 


P1  =  2P  -  P'1  +  ^  [(P1  -  P0)C^  -  <P0  -  P  3 ) ] 

hl 


+  ^-t(P2  -  P0)O2  -  <P0  -  P4>“4] 


+  r1-  [B?s  -  623s2] 

h ,  1  z  i  z 


+  [$2s  -  ] 

h„  4  x  2  x 


(4.2.19) 


1  -1  At' 

S  =  2S  -  S  +  — 


t(Sl  -  V61  -  <S0  -  S3)B3] 


[<S2  -  V62  -  (S0  -  S4>@4] 


At 


[Vx  -  VJ  +  hT  [Vz  -  xipz]<4-2-20) 
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We  note  that  the  boundary  terms 


h“  (f3]S 
h ,  1  z 


63sz) '  ir  <64s  -  b^s  >  , 

3  z  h^  4  x  2  x 


(X  P  -  X  P  )  ,  —  ( X  P  -  X  P  ) 
2  x  4  x  h  3  z  lz 


(4.2.21) 


represent  the  conversion  terms  SV  to  P  and  P  to  SV. 

The  finite  development  of  those  terms  have  the  form 


( efs  -  e^S  )  =  — --r-  ($2  -  3?)  (s  -  S  )  (4.2.22) 

h  lz  3z  2hh  1  3  2  4 

I  J-  £ 


12  2  1 
(B^S  -  (Ts  )  % 
h  _  4  x  2  x 


2h  h  (B42  *  e22)<Sl  -  S3> 


1  2 


(4.2.23) 


4.2.b  Development  by  Taylor  Expansion 


Equation  (4.2.2)  and  (4.2.3)  yield 


r  2^2  -I  r  2  2,^  „ 

F.D.[a  V  p]  -  [a1  -  Ot  3  ]  Dxp 


2  2  .  'U 

-  (a_  -  a. )  D  P 

2  4  z 


=  FD  (t5 ) 


(4.2.24) 


fd[62V2s]  -  <B2 


2 

33)dxS 


-  ( 3?  -  3^)D  s  =  F.D(S) 

2  4  Z 


(4.2.25) 
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Since  the  conversion  terms  are  identical  to  expressions 
(4.2.15)  to  (4.2.18)  the  development  leads  to: 


F  D { V  2 }  E  FV {  (  V)  ? } 


The  development  of  equations  (4.2.25)  and  (4.2.26)  is  iden 
tical  to  equations  (4.2.19)  and  (4.2.20)  obtained  by  varia¬ 
tional  principle.  i.e.  After  taking  into  account  the 
boundary  conditions  (4.2.11)  to  (4.2.14),  (4.2.25)  and 

(4.2.26)  yield 


F D ( P )  =  F  D [a2 V2P]  +  $2S  -  R2s 

x  z  z  x 


(4.2.26) 


FD(S)  =  F  DfB2V2S]  +  X  P  -  X  P 

fz  X  -X  z 


(4.2.27) 


which  is  identical  to  (4.2.19)  and  (4.2.20)  or  in  the  form 


r\j 


'U 

S 


=  a2^2P  + 

2^ 
a  P 

X  X 

2^ 

+  a  P 

Z  z 

+  e2s 

X  z 

r,2^ 

-  3  S 
z 

=  $2V2s  + 

3  s 

X  X 

2^ 

+  3  s 

z  z 

+  X  P 
^rz  x 

-  X  P 

rx 

(4.2.28) 


(4.2.29) 


where  ^  represents  the  homogeneous  development  of  section 
2.4.  We  have  expressed  the  form  (4.2.28)  and  (4.2.29)  not 
in  a  computational  form  but  for  qualitative  purposes. 
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4.2.c  Reflection  and  Conversion  Terms  Considered  as  Source 


We  notice  that  in  equation  (4.2.28)  for  instance 
2  2 

axpx'  azp2  represents  the  reflection  and  refraction  terms. 
2  2 

(£>xsz  -  f^S^)  represents  the  SV-*P  conversion  term 


a.  'u 

(X  p  “  _A  P  )  represents  the  P-^SV 


^  Z  X  ^  X  z 


conversion  terms. 

From  (4.2.28)  and  (4.2.29)  we  note  the  symmetry  of 
reflection  terms  and  conversion  terms. 


Be  sides , 


equation  (4.2.29)  shows  that  even  in  a 


homogeneous  fluid  SV  may  be  generated  but  does  not 


propagate;  since  (4.2.27)  can  be  written 


S  =  curl  K 


(4.2.30) 


with 


curl  K  =  X  P  -  X  P 
z  x  x  z 


(4.2.31) 


But 


S  ( M )  E  0  V  F(M)  =  0 

V  M  £  ft  c  1R2  x  [ot] 


In  (4.2.30)  we  have  considered  the  transformation  of  a 
source,  i.e.  as  a  body  force.  Then  in  equations  (4.2.26)  , 
(4.2.27)  ,  the  reflection  and  transmission  terms  can  be 
considered  as  sources  (or  body  forces)  due  to  a  discontinuity 


(parameter  discontinuity) ,  those  body  forces,  satisfy 
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of  course  the  B.C.  Then  from  (4.2.26)  and  (4.2.27)  we  can 
write  the  reflection  and  conversion  terms  considered  as 
body  forces 


V.K  =  a2P  +  a2P  +  $2S  -  $2s 

XX  zz  xz  zx 


(4.2.32) 


Vak  =  B2s  +  $2s  +  A  p  -  A  P 

XX  Z  Z  ^~Z  X  ^rx  z 


(4.2.33) 


By  analogy  with  the  equation  of  displacement  motion  we 
can  write  Kx ,  Kz  as  the  components  of  the  reflection 

and  conversion  coefficients 


2  2  2 
Kx  =  a  u  +  A  w  +  3  u  +  3  w 

X  X  x  Z  Z  Z  Z  X 

e 


(4.2.33) 


Kz  =  a2w  +  A  u  +  3  w  +  3  u 

zz  ^-zx  XX  xz 


(4.2.34) 


But 


curl  K  =  Kx  -  Kz 
z  x 


Expressing  curl  K  with  the  help  of  (4.2.33)  and  (4.2.34)  , 
and  taking  into  account  that 


P  =  U  +  W 
x  z 


s  =  u  -  w 

z  x 


* 
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as  well  as  the  B.C 


+ 

a  =  o 

xz 


we  obtain 


2  2 

curl  k=3s  +  3  s  +  A  p  -  A  p 

XX  zz  zx  xz 

which  is  the  relation  (4.2.33)  ,  which  can  be  considered 

as  a  verification  of  the  assertions  given  in  this  section. 
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4.3  "HETEROGENEOUS  MEDIA  SOLUTION"  OF  EQUATIONS  OF 
COMP  RE  S  S I ONAL  AND  ROTATIONAL  DEFORMATION 

4. 3. a  Development  of  the  so  called  "heterogeneous  media 
solution" 

Following  Landers  and  Claerbout  (1972)  ,  Grant  and 
West  (1965)  ,  we  consider  the  law  of  motion  for  the  dis¬ 
placement  vector  U,  i.e. 


pu  =  (A+  y)V(V.u)  +  VAV.u 


-  pVaV  u  +  Vy  a (Vau) 


+  2  ( VyV ) u 


(4.3.1) 


(Karal  and  Keller,  1959) 


Letting 


V.U  = 


P 


Vau  =  s 


We  obtain  after  taking  the  div.  and  curl  of  (4.3.1) 
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P  =  a2  [p  +  p 

XX 


zz 


2  2 

)  +  2a  P  +  2a  P 

XX  z  z 


(4.3.2) 


s  =  $2  [S  +  S 


XX 


z  z 


)  +  2  $  S  +  2^  S 


XX  z  z 


(4.3.3) 


(Landers  and  Claerbout,  1972) 


4 . 3 . b  Remarks 

We  see  that  if  the  "heterogeneous  media"  method  is 
correct  we  shall  have  two  important  properties 

i)  From  (4.3.2)  and  (4.3.3)  the  conversion  terms  are 
independent  of  A,  i.e.  waves  can  be  converted  only  if  there 
is  a  discontinuity  in  shear  velocity. 

ii)  If  we  consider  a  fluid  equation  (4.3.2)  becomes 


instead  of 


2 

+  a  P 
z  z 
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i.e.,  the  reflection  coefficients  will  be  twice  the  reflec¬ 


tion  co 


of  a  scalar  wave  equation  or  a  sound  wave 


Therefore  the  system  (4.3.2),  (4.3.3.)  is  obviously  not 


correct . 


4.4  General  of  P-*SV,  SV->-P  converted  waves 

Elastic  wave  motion  is  governed  by  equations  (4.2.29), 
(4.2.3),  i.e. 

-  variation  in  X  generate  P^-SV 

-  variation  in  y  generate  SV-^P  . 

Figures  (4.1a)  to  (4. Id)  represent  snap  shots  of  P  S  V 

converted  wave  (P  source)  for  a  simple  model,  (Fig.  4.1); 

the  source  having  a  Gaussian  dependence  in  a  domain  of 

2 

definition  £2  =  P  ±  2h  c  R  x  [0,t)  and  where  the  spatial 

s 

dependence  S  (  X  )  of  the  source  trends  smoothly  to  zero 
at  the  boundary  of  its  domain  of  definition. 

The  model  consists  of  two  layers  of  7000'/sec  and 
lOjOOO’/sec.  The  wave  field  computed  with  a  coarse  grid 
is  represented  for  qualitative  purpose  and  shows  the 
relative  importance  of  converted  shear  wave. 


Fig.  4.1c  represents  the  shear  wave  field,  and 
comparison  with  the  P-*SV  reflection  coefficient  curve  as 
a  function  of  the  angle  of  incidence  shows  the  corre¬ 
lation  between  the  theoretical  curve  and  the  synthetic 


way e  field. 


. 


REFLEXION  AND  REFRACTION 
OF  A  COMPRESSIONAL  WAVE 
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o 

o  W 

LU  CO 


ro 


Fig.  4.1.  Reflection  and  conversion  of  a  compr e s s i ona 1  wave 
4.1. a  One  layer  model 
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yooo'/se^ 


lOOOO/s  l 


SNAP  SHOT;  WITH  P  SOURCE 


P  WAVE  ONLY 


( 1~  —  -  -90  xj ec.J 


P  DIRECT 


P  REFLECTE 


P  TRANSMIT 


Fig.  4.1.b  Wave  field  of  reflected  and 
transmitted  P  wave 
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SNAP;  SHOT;  WITH;  P:  SOURCE 
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s’-’Ai'  K!i° T:  WITH'  P.  S O U R C E 

*  •  9o  S£^.)  " "  ~  ' 


-  -  3  So  - 


-  7000^ec 


10  000%c 


=VT 


-P  DIRECT 

f 

j 

P  REFLECTED  j 

* 

SV  REFLECTED ; 


SV  l  RARSMITED 


P  TRANSMITED 


Fig .  4 . 1 . d  Total 


vave  field. 


CHAPTER  V 


Cone lus ion 


The  main  points  of  this  research  which  are  original 
can  be  summarized  as  follows: 

i)  The  distribution  of  the  parameters  at  the  boundary. 
Since  we  can  consider  a  measure  as  the  inverse  of  a  partial 

(0^  ) ,  it  is  normal  to  consider  the  parameters 

as  specific  to  the  operator  domain  (in  this  case  the  domain 
of  the  second  derivative) .  An  important  consequence  being 
that  the  parameters  can  vary  with  the  direction  and  diagonal 
boundary  can  be  used. 

ii)  The  difficulty  of  introducing  the  boundary  condi¬ 


tions  inside  the 


equation  for  the  direct  method 


has  been  solved  simply  by  adapting  the  Taylor  expansion  to 
the  specificity  of  the  parameters.  We  note  that  if  the 
boundary  conditions  do  not  vanish  when  introduced  in  the 
differential  equation,  the  so  called  "heterogeneous  media" 
method  (Kelly  et  al,  1976)  cannot  be  used,  which  is  the 
case  of  most  problems  except  in  the  SK  or  the  scalar  wave 
equation.  The  direct  scheme  is  valid  for  any  case  of 
boundary  conditions  and  can  be  considered  as  a  direct 
generalised  scheme  for  heterogeneous  media. 
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iii)  The  Dirichlet  integral  and  the  energy  method  give 
us  a  classical  way  of  introducing  the  boundary  conditions. 
The  z  transform  and  the  choice  of  a  trial  function  which 
has  the  same  norm  as  for  the  direct  method  (in  Sobolev 
space  which  is  a  Hilbert  space)  allow  us  to  obtain  a 
finite  development  which  is  comparable  to  the  direct  form. 
The  identity  of  the  schemes  and  their  uniqueness  in  the 
two  methods  developed  shows  the  problem  is  well  posed. 

iv)  The  P  SV ,  and  SV  P  transformations  have  been 
studied  through  the  divergence  and  rotation  or  curl  of  the 
displacement  vector  equation.  The  direct  and  the  finite 
variational  developments  have  been  successively  applied 
and  lead  to  identical  schemes.  Those  solutions  correlate 
with  the  displacement  vector  solutions  seen  above.  The 
parameter  discontinuity  is  considered  as  a  source  of  the 
scattered  and  converted  fields  in  the  displacement  vector 
equation.  The  divergence  and  curl  of  those  forces  give  the 
scattered  and  converted  field  in  terms  of  the  divergence  and 
curl  of  the  vector  displacement  equations.  In  this  case  als 
the  comparison  with  the  so  called  "heterogeneous  media" 
method,  gives  the  same  conclusion  as  before.  Although  the 
results  are  not  concordant  to  the  one  obtained  by  Landers  et 
al  (1972) ,  and  Grant  and  West  (1965) ,  the  evidence  led  us 


142 


to  believe  that  the  SV-*P  and  P->SV  conversions  are  given 
respectively  by  discontinuity  of  the  parameters  y  and  A 
(and  not  only  y  as  stated) .  This  result  is  important: 

-  as  an  interpretative  tool  in  seismic  exploration. 

-  to  investigate  the  bright  spot  problem  (Stoffa  et 

al  ,  1976)  . 

v)  The  secondary  results  of  this  research  are: 

-  the  representations  of  the  source  as  the  causal 
part  of  the  development  of  a  force.  This  underlines  the 
reciprocal  properties  of  the  wave  equation  and  is  very 
useful  in  the  numerical  solution. 

-  The  finite  difference  development  of  the  SH  or 

scalar  wave  equations  at  normal  incidence  for  V _ t  =  1 

h 

is  an  exact  solution.  That  is,  the  simple  finite  difference 
development  is  identical  to  the  analytical  plane  wave  solu¬ 


tion  or  the  matrix  method. 


■ 
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